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Preface 


Since 2014, we have organized an annual summer school in computational physiol- 
ogy. The school starts in June each year and the graduate students spend two weeks 
in Oslo learning the principles underlying mathematical models commonly used in 
studying the heart and the brain. At the end of their stay in Oslo, the students are 
assigned a research project to work on over the summer. In August the students travel 
to the University of California, San Diego to present their findings. Each year, we 
have been duly impressed by the students’ progress and we have often seen that the 
results contain the rudiments of a scientific paper. 

Starting in the 2021 edition of the summer school, we have taken the course one 
step further and aim to conclude every project with a scientific report that passes 
rigorous peer review as a publication in this new series called Simula SpringerBriefs 
on Computing — reports on computational physiology. 

One advantage of this course adjustment is that we have the opportunity to intro- 
duce students to scientific writing. To ensure the students get the best introduction 
in the shortest amount of time, we have commissioned a professional introduction to 
science writing by Nature. The students participate in a 2-day Nature Masterclasses 
workshop, led by two editors from Nature journals, in order to strengthen skills in 
high quality scientific writing and publishing. The workshop is tailored to publica- 
tions in the field of computational physiology and allows students to gather real-time 
feedback on their reports. 

We would like to emphasise that the contributions in this series are brief re- 
ports based on the intensive research projects assigned during the summer school. 
Each report addresses a specific problem of importance in physiology and presents 
a succinct summary of the findings (8-15 pages). We do not require that results 
represent new scientific results; rather, they can reproduce or supplement earlier 
computational studies or experimental findings. The physiological question under 
consideration should be clearly formulated, the mathematical models should be de- 
fined in a manner readable by others at the same level of expertise, and the software 
used should, if possible, be made publicly available. All reports in this series are 
subjected to peer-review by the other students and supervisors in the program. 
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A 
Chapter 1 ‘patos 
An Automated Cardiac Constitutive Modelling 
Framework with Evolutionary Strain Energy 
Functions 


Kristin Ludwicki*, Leto L Riebel*, Sophia Ohnemus*, Frida M E Westby, 
Nickolas Forsch, and Gabriel Balaban 
* These authors contributed equally to this work 


Abstract Heart disease is the leading cause of mortality worldwide. Many cardiac 
diseases are associated with altered elastic energy relations of the heart tissue. How- 
ever, the strain energy functions describing these characteristics are limited since 
they were designed manually for highly specific experimental setups. In this study, 
we develop CHESRA (Cardiac Hyperelastic Evolutionary Symbolic Regression Al- 
gorithm), an automated constitutive modelling framework, to derive cardiac elastic 
strain energy functions directly from experimental data. Our results indicate that 
CHESRA finds functions that reproduce mechanical tissue properties from exper- 
imental data whilst controlling function complexity. Our novel approach has the 
potential to find strain-energy functions that fit to various experimental data sets 
and may contribute to automatically building mathematical models to understand 
clinical observations of heart diseases. 
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2 An Automated Cardiac Constitutive Modelling Framework 


1.1 Introduction 


Cardiovascular disease is the leading cause of mortality worldwide, accounting for 
around 19 millions deaths in 2020 [1]. In many cardiac diseases, the mechanical 
and structural properties of the heart are altered. Examples include hypertrophic 
cardiomyopathy, which is characterized by a thickening of the heart walls [2] and di- 
lated cardiomyopathy, a disease which causes the ventricular walls to thin and stretch 
[3]. Another example is myocardial infarction, commonly known as a heart attack, 
during which oxygen and nutrient supply to the cardiac tissue is interrupted [4]. 
Myocardial infarction results in altered mechanical loading conditions [5], structural 
changes, increased tissue stiffness, and ultimately ventricular dysfunction [6]. A bet- 
ter understanding of these pathological mechanisms may improve clinical diagnosis 
and treatment options. 

To characterize the elastic properties of the heart, stress-strain relations of the 
myocardium have been extensively studied experimentally. Assessment of passive 
shear properties of pig ventricular myocardium by Dokos et al. [7] indicates that 
the myocardium presents as an orthotropic material since the shear response differs 
along three mutually orthogonal directions. Demer and Yin [8] and Yinet al. [9] mea- 
sured the stress-strain relationship of passive non-contracting canine myocardium for 
simultaneous biaxial stretching and observed highly nonlinear, anisotropic, and vis- 
coelastic behavior. Novak et al. [10] performed more biaxial experiments in canine 
myocardium from different regions of the left ventricle and found that qualitatively 
the mechanical properties were similar, but quantitatively they differed. For example, 
samples from sub-endocardium and sub-epicardium of the left ventricular free wall 
tended to be more rigid than samples from the mid-myocardium. 

In order to mathematically represent the mechanical properties of the heart tissue, 
a number of strain-energy functions (SEF) have been proposed. These SEF describe 
the potential energy density of the heart tissue depending on its deformation and 
have been developed based on a few key considerations. Considerations include 
structural metrics and other experimental results, as well as assumptions on physical 
properties of the system. Although the myocardium appears to be a viscoelastic 
material, the relaxation time of the viscoelastic response is long compared to the 
cardiac cycle [11]. Therefore, most SEF, with exceptions as presented in [12], model 
it as hyperelastic for simplicity reasons. Examples for SEF based on the assumption 
of transverse isotropy are given in [13, 14], whereas [11, 15, 16, 17] introduce 
orthotropic models. However, creating these SEF is labor intensive as they must 
be derived manually and, thus, only a limited amount of experimental data can be 
considered. 

Evolutionary algorithms have the potential to automatically find functions that 
describe experimental data with minimal human guidance. This has been applied 
already to various research fields; for example, Doglioni et al. used evolutionary 
polynomial regression (EPR) to build a function for the relationship between air 
and water temperature [18]. Their method not only considered suitability of the 
discovered function to describe the given data, but also aimed to minimise the 
complexity of the formula. Montes et al. used a modified EPR framework to create 


1 An Automated Cardiac Constitutive Modelling Framework 3 


self-cleansing models for new sewer systems [19] and several different machine 
learning methods, including EPR, were used by Jamei et al. to model the density 
of hybrid nanofluids [20]. Furthermore, Javadi et al. established an EPR method to 
build constitutive equations describing the deformation of bridges and tunnels [21]. 

In this work we present the development of a Cardiac Hyperelastic Evolutionary 
Symbolic Regression Algorithm (CHESRA) that derives hyperelastic SEF for ven- 
tricular myocardium. In contrast to EPR algorithms, CHESRA is not restricted to 
polynomials, but searches the whole space of mathematical functions for SEF de- 
scribing the given experimental data. Our results highlight the potential to automat- 
ically find SEF which reproduce experimental findings of cardiac elastic properties 
of different species. 


1.2 Methods 


In the following subsections we outline the individual components of CHESRA. 
First, in 1.2.1 we briefly introduce the physical principles of hyperelastic continuum 
mechanics that we use to define the SEF. Then, in 1.2.2 we outline how we represent 
the SEF computationally as function trees. Subsequently, we explain how we fit the 
SEF to experimental data in 1.2.3. Lastly, we go through the individual steps of the 
evolutionary symbolic regression (ESR) algorithm in detail in 1.2.4. 


1.2.1 Basic Principles of Hyperelastic Continuum Mechanics 


In order to define our SEF, we build upon the framework developed by Holzapfel 
and Ogden for cardiac hyperelasticity [11]. Here, the structure of the myocardium 
is characterized by three orthonormal basis vectors: the fibre axis f associated with 
the prevailing cell orientation, the sheet axis s defined as the direction in the plane 
of the cell sheets perpendicular to the fibre direction, and the sheet-normal axis n 
perpendicular to f and s. 

Any point within the material can be described by a vector X = ))j-+,5  Ci€i; 
where e; denote the basis vectors of the f, s, and n axis. Upon deformation, the 
position of each point changes and the new position can be described by a vector 
x. The fundamental quantity describing such a deformation is the deformation gra- 
dient F, defined by F;; = 0x;/0X; with i, j = f,s,n and where J = det F = 1 for an 
incompressible material such as the myocardium. The right Cauchy-Green tensor 


C=F’F, (1.1) 


and the Green-Lagrange strain tensor 


1 
E=5(C-D, (1.2) 
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are associated with F, with I being the identity tensor. 

Since these quantities do not only account for deformation but also for rotation 
and translation, in the following, we consider the invariants of C which are defined 
as 


1 
I, =trC, haz [77 -trC] b = detC, 
14; =e} (Cei), Is; =e; (C’e;), Igi; =e; (Ce;), (1.3) 


with i + j € {f,s,n}. Defining a SEF based on these invariants ensures objectivity. 
The SEF proposed by Holzapfel and Ogden [11] is given by 


w=, PLC -3)] p2 py Plb: Ca- 171-1} 
L lexplbysiy.1— 1}. (1.4) 


+ 
2b fs 


Here, a, b, af, bf, as, bs, afs, and bys are material parameters chosen to either fit 
the corresponding Cauchy stress tensor 


ror FS ——_— (1.5) 


to the shear data from Dokos et al. [7], or the associated second Piola-Kirchhoff 
stress tensor 


S=JF ‘oF? (1.6) 


to the biaxial stretch data by Yin et al. [9] (see Section 1.2.3.1 for details about these 
experimental data sets). 


1.2.2 Representation of the Strain-Energy Functions as Function Trees 


In CHESRA, SEF are represented as a list of nodes composing function trees. Details 
on function tree implementation and use cases can be drawn from [22, 23]. In our 
implementation, each node in the function tree has a type, a value, up to one parent, 
and up to two children (one left and one right). Only the root node of the function 
tree has no parent, further explained in 1.2.4.1. Nodes with no children are called 
*leaves’. Sub-trees within the function tree represent individual terms of the SEF. 
Nodes are either a symbol, an operand, or a pre-operand. Nodes of type symbol 
can either be a constant material parameter or an invariant. Material parameters are 
represented by the symbols {p1, P2, P3, P4, P5, P6.P7,Ps}- These material param- 
eters are later replaced by a constant real number to fit the SEF to experimental 
data, as described in 1.2.3. Similarly, invariant nodes are assigned a symbol from 
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the set of invariants {1), Jy, 13, Tap, 14s, Lan, 15 ¢,155,15ns1g¢s,18fn>Igns}, as defined in 
Equation 1.3. Operand nodes can either have the value + (addition), — (subtraction), 
+ (division) or * (multiplication); and pre-operands are either pow (to the power of 
2), exp (euler’s number ’e’ to the power of) or - (negation). While operands require 
two terms, for example x + y, a pre-operand requires only one, for instance e*, which 
is written before the term in many coding languages, such as in python: pow(x, 2). 
An example function tree is shown in Figure 1.1. 


Ų = pı * (2 + lf) 


T Operand 
Constant Pre-Operand Variable +4 


Fig. 1.1: Example SEF y and corresponding function tree with node types indicated. 
The function tree is converted to a real function y by applying an in-order traversal. 


If a node of type symbol, i.e. either a material parameter or an invariant, has a 
parent node, it has to be either an operand or a pre-operand to be a valid function. 
Similarly, symbol nodes can only have children of type operand. Pre-operands are 
applied to terms, and hence can only have exactly one child, which has to be of 
type symbol or pre-operand. Their parent has to be either an operand or another 
pre-operand. Operands have a parent and exactly one child, both of which cannot be 
another operand. If an operand is a right child, its child needs to also be a right child 
and vice versa, to avoid two symbols next to each other. These rules are summarised 
in Table 1.1. We apply brackets only around the terms pre-operators are applied to, 
all other operations are carried out according to their priority, e.g. + before +. 

Using the rules summarized in Table 1.1, to convert our function trees into 
functions we apply an in-order traversal, where any node is read out as: value of its 
left child, its own value, value of its right child. Figure 1.1 shows an example of 
a function and one of its possible corresponding function trees, depending on the 
order the terms are created in. Note that since our algorithm does not perform any 
searches within the tree, there is no need for it to be balanced. 

Representing SEF as function trees has several benefits. Firstly, individual terms 
can be identified as sub-trees, allowing for pre-operands to be applied to parts of the 
equation only. This would be more difficult, for example, in a linked list structure. 
Secondly, terms can be manipulated, swapped, and deleted easily during mutation 
and mating, which are explained in 1.2.4.2. 
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Type Value Rules 


Symbol Invariant (variable), -Parent has to be an operand or pre-operand 


material parameter (const) | -Child has to be an operand 


Operand +, —, *, + -Has exactly one parent of type symbol 


-Has exactly one child of type symbol/pre-operand 


Pre-Operand | e*, x^, —x -Parent has to be an operand or pre-operand 


-Child has to be a symbol or pre-operand 


Table 1.1: Summary of node types, possible node values, and node connection rules 
within the function tree. 


1.2.3 Fitting the Strain-Energy Functions to Experimental Data 


For each SEF y we fit the material parameters pn to experimental data sets in a 
separate step, in order to ensure that only the material parameters differ for different 
experimental setups while the general form of the SEF is the same. In this work we 
consider the shear data of Dokos et al. [7] and the biaxial stretch data by Yin et al. 
[9], which we briefly describe in the following. 


1.2.3.1 Experimental Data under Consideration 


Dokos et al. [7] measured the shear stress, corresponding to the Cauchy stress tensor 
components o;; versus the amount of shear y in a cube of pig left ventricular 
myocardium sheared in the fs, fn, or ns plane (see Figure 1.2a). They considered 
different shear modes (ij) referring to shear in the ij plane in j direction with 
i#je{f,s,n}. 

The biaxial stretch data from Yin et al. [9] was collected in canine left ventricular 
tissue which was stretched simultaneously in the fibre and sheets direction. Figure 
1.2b shows the measured second Piola-Kirchhoff stress Sz versus strain Eff and 
respectively Sss versus Ess for three different constant strain ratios r = Ey ¢/Ess. 

For both data sets, we used the digitizer software WebPlotDigitizer [24] to digitize 
the data from the graphs shown in [7] and [9]. 


1.2.3.2 Fitting to the Experimental Data 
In order to fit the material parameters pn of a SEF w to the shear data set we 
calculated the predicted Cauchy stress tensor components oj; according to Equation 


1.5. Then, we used the Python package Imfit [25] to minimize the residuals, 


resk = Sij (Vk) — Sij (Yk), (1.7) 
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(a) Shear data (b) Biaxial stretch data 


Fig. 1.2: Experimental data sets under consideration. (a) Plot of the digitized data 
collected by Dokos et al. [7] in pig left ventricular myocardium for shear stress 
oj; versus the amount of shear y for shear (ij) in the ij plane in j direction with 
i+ j €{f,s,n}. Results for (nf) and (ns) shear are identical. (b) Visualization of 
the data by Yin et al. [9] for biaxial loading in the fs plane of canine left ventricular 
myocardium. The left plot shows the stress S$ p f versus strain E + ¢ in the fibre direction 
and the right plot S,, versus Ess in the sheets direction for three different strain ratios 
r = Eşff/Ess =2.05 (triangles), 1.02 (squares) and 0.48 (circles). In both cases we 
used WebPlotDigitizer [24] to digitize the data shown in [7] and [9]. 


between experimental data points s;; and the corresponding calculated stress tensor 
component o;; for all shear modes (ij) and all amounts of shear yx in the data set 
simultaneously. 

Analogously, for fitting the material parameters to the biaxial data set we cal- 
culated the diagonal components of the second Piola-Kirchhoff stress tensor S;; 
according to Equation 1.6. We used the same Python package for simultaneously 
minimizing the residuals 


resk = Sii( Eiji.) — Si ii); (1.8) 


with i = f,s, for all Green-Lagrange strain tensor components Eji k and strain ratios 
r = Eff/Ess in the data set. 


1.2.4 Evolutionary Symbolic Regression Algorithm 


To develop CHESRA, we build upon the traditional ESR framework. First, popula- 
tions of random SEF are created (see Section 1.2.4.1). Then, every equation within 
each iteration, or generation, is assigned an error-score based on a user-designed 
fitness function. Our fitness function assesses the complexity of the SEF and how 
well it fits to the shear and/or biaxial data set (see Section 1.2.4.3). The individuals 
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with the best fitness are then selected, mated to generate a new population, and mu- 
tated to ensure diversity for the next generation (see Section 1.2.4.2). This multi-step 
process continues until the maximum number of generations is reached. A schematic 
representation of CHESRA is illustrated in Figure 1.3. 


INITIALIZATION | aE a oe 
Generate a popula- + í paaa `- START CHEPRA 
tion of function trees <=. a 


t 


build function trees 


+ 
(EVALUATE FITNESS 
Select training data and up- 
date fitness function based on | — —_ 
the number of training data 
sets and function complexity 


No 


_—_———_£ A —— 


/ load in training ~~ /~ count equation ~N 
( data set(s) and { length and number 
N fit functions / of pre-operators 


calcula! 3 $ aji 
( calculate 


__ for each fur complexity for 
each function 


calculate total 
fitness (SSE + 

Complexity) for 
cach function / 


ł 


SELECTION 
Sets of two equations are ran- 
domly chosen and the individual > aming yot are rorideniy selected 10 m_e 


selected for mutation at a generation? / 
mate at a probability set by the user. á 


with the lower error score is probability set by the user. 


then added to the mating pool. 


MATING ) MUTATION yo 
Sets of two equations from the Equations can be randomly Last 


Fig. 1.3: The main workflow of CHESRA. 


1.2.4.1 Setup and Initialisation 


To initialise SEF and their corresponding function trees, a maximum function length, 
set of invariants, and material parameters are defined. First, our algorithm creates 
a root node, which is always an invariant or a material parameter. Second, it picks 
a random length within the given maximum and extends the function tree until this 
chosen length is reached. The function can either be extended by adding a pre-operand 
in front of a term, or by attaching an operand and a new symbol (i.e., an invariant 
or material parameter) to an existing symbol. For the latter, the selected existing 
symbol may have at most one child. Details concerning function development can 
be found in Section 1.2.2. 
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1.2.4.2 Mutation and Mating of Strain-Energy Functions 


We mutate SEF by randomly selecting a node from the corresponding function tree 
to be mutated. The selected node is then replaced by a randomly chosen element of 
the same type (see 1.4). 

To mate two parent SEF and create two new children, our algorithm creates a 
deep copy of each parent function tree, chooses one random sub-tree per copy, and 
then swaps them over. To make sure the resulting children SEF are valid functions, 
we restrict the chosen sub-trees to be proper terms and hence start with a symbol 
or pre-operand. Alternatively, one may also choose to ensure the sub-trees to be 
swapped to start with a node of the same type. Our chosen sub-trees may be of 
different size and hence children SEF may be of different lengths than their parents 
(see 1.5). 

Mutation and mating are important to increase population diversity in CHESRA. 
However, too much diversity could also prohibit convergence towards an optimal 
solution. Therefore, we set rates that will dictate the probability of a given individual 
to mate, Pmate, and/or mutate, Pmutate. These rates were set tO Pmate = 20% and 
Pmutate = 80% in all CHESRA experiments. 


Y= pr * (Pot lay)? Y =p * (p2— lap)? 


Fig. 1.4: An examplary function tree representing a SEF y is mutated by choosing a 
random item and replacing it with a random value of the same type, as highlighted 
by the orange boxes. 


1.2.4.3 Scoring Fitness of Strain-Energy Functions 


We designed a fitness function to ensure that CHESRA can penalize against poor 
fitting SEF with high complexity. Therefore, each SEF, or individual, in each gener- 
ation is evaluated based on its complexity and how well it reproduces the results of 
a specific experimental data set. It is this evaluation that allows the best individuals 
to be selected for mating, mutation, and population of the following generation. 
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Fig. 1.5: Mating between two parent SEF to generate two children SEF is imple- 
mented by swapping two random subtrees from the corresponding function trees as 
highlighted in the orange boxes. 


To evaluate the fitness of each equation, first we quantify its complexity as the sum 
of the function length Leq (i.e., the number of nodes it consists of), and the number of 
pre-operators pre-op that it contains. The pre-operator count allows for penalization 
against nestedness, or the creation of functions with pre-operators within themselves. 


Complexity = leq +Npre-op- (1.9) 


Second, we assess how well each individual SEF reproduces experimental data. 
For this, the material parameters of each equation are first fitted to the experimental 
data as described in 1.2.3.2. Then, for each equation the sum of squared errors (SSE) 
is calculated as the sum of the squared residuals res, between the experimental data 
points and the fit, 


SSE =) res, (1.10) 
k 


with res, as defined in Equation 1.7 for the shear data or in Equation 1.8 for the 
biaxial stretch data. 
The final fitness score is then given by 


Fitness = (œ x Complexity) +SSE, (1.11) 


where the hyperparameter œ ensures appropriate balance between fitting and com- 
plexity. Numerical experiments conducted to find a suitable œ can be found in Section 
131; 

Where CHESRA is given N data sets, the fitness function can be extended to 


Fitness = (a* Complexity) + X SEN. (1.12) 
N 


1 An Automated Cardiac Constitutive Modelling Framework 11 


Hence, in the following section 1.3.1 we used Equation 1.11 in order to find a 
SEF that can reproduce the shear data alone, while in section 1.3.2 we used Equation 
1.12 to consider both the shear and the biaxial data set. 


1.3 Results 
1.3.1 Trial 1: SEF Describing the Shear Data 


In our first trial, we ran CHESRA while fitting the SEF to the shear data set from 
Dokos et al. [7]. We chose the hyperparameter a of the fitness function (see Equation 
1.10) by running CHESRA with a = 1075, 1073, 107!, 1, 2, and 5 for 50 generations 
and 100 individuals each. 1.6 indicates that a value of a =0.1 lead to the best result as 
it allowed for SEF with low SSE and minimal complexity. Using this hyperparameter, 
the fitness score of the best SEF gradually decreased over the generations until a 
plateau was reached (see 1.7). 

In this run, the SEF with the lowest fitness score derived using CHESRA repro- 
duces the shear data well (see 1.8a) and is given by 


Wy = exp (pi Iss + prts¢ +2h — Ign — p3). (1.13) 


In order to assess whether CHESRA is reproducible, we repeated this trial twice, 
with the best SEF given by 


y2 = —-CUs ¢ — pı} + pals p(s + ps)? + las), (1.14) 
y3 =exp(pilsp t+ (pot p3—l)Igps + palag +l) +212 — Ian + ps). (1.15) 


The equations derived from the three trial runs have unique material parameters 
that can be found in 1.2. 1.8 shows a good fit to the shear data for all three SEF found. 


Trial run | pı p2 P3 P4 P5 


1 0.77 | 2.00 | 9.92 | - - 
2 1.47 | 0.20 | 0.50 | - - 
3 0.86 | 1.44 | 0.04 | 1.38 | 510.06 


Table 1.2: Material parameters to fit the SEF defined in Equation 1.13 (trial run 1), 
1.14 (trial run 2), 1.15 (trial run 3) to the shear data set by Dokos et al. [7] 
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Fig. 1.6: SSE vs. complexity for different hyperparameters a in the fitness function 
(see Equation 1.11). Six runs of CHESRA were completed in which each was run 
with one of the unique œ values listed in the legend. For experiment 1, we chose 
a =0.1 since it lead to SEF with low SSE and function complexity. 
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Fig. 1.7: Evolution of 100 individual SEF for 50 generations. Each point represents 
the fitness of the best SEF from each generation. 
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Fig. 1.8: In our numerical experiment 1, CHESRA found a SEF that reproduces the 
shear data of Dokos et al. [7] in all 3 trial runs. Markers represent the experimental 
data points whereas solid lines show the fit of the SEF derived using CHESRA. 


1.3.2 Trial 2: SEF Describing the Shear and Biaxial Data 


As a next step, we extended CHESRA to find a SEF that reproduces the results of 
multiple experimental data sets. Therefore, we ran CHESRA with the same setup as 
in trial 1, but updated the fitness function to incorporate the SSE of both the shear 
and the biaxial data set, as described in Equation 1.12. Moreover, in this numerical 
experiment we used a = 0 in order to assess quality of fit without penalizing for 
equation complexity. In addition, a = 0 allowed convergence in a shorter number of 
generations as fitting to two different data sets is a challenging optimization problem. 
The SEF found using CHESRA is given by 


1 
w= |piprlap— Tp + = - las [P1 Pzlsnla4f — pa p4h + ps + pot Pr 
5n 


2 
- paps pi ps py ; (1.16) 


where the material parameter values are provided in Table 1.3. Figure 1.9 shows that 
the equation derived in experiment 2 has a better fit to the biaxial data than to the 
shear data. 


Data | pı P2 P3 | P4 | Ps | Po | P7 


Shear | 0.00 | 433.64 | 33.52 | 0.31 | 1.85 | 3.08 | 4.95 
Biaxial | 0.49 | 1.34 1.24 | 1.11 | 2.56 | 0.04 | 2.06 


Table 1.3: Material parameters for fitting the SEF defined in Equation 1.16 to the 
shear and biaxial data sets. 
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Fig. 1.9: In experiment 2, CHESRA was extended to find a single SEF aimed to 
reproduce both the shear data from Dokos et al. (a) and biaxial stretch data from 
Yin et al. (b). Solid lines show the fit of the SEF found using CHESRA to the 
experimental data represented by the markers. 


1.4 Discussion 


CHESRA automatically generates functions which reproduce experimental data of 
cardiac hyperelastic properties. The results of trial 1 indicate that CHESRA can 
derive SEF that fit the shear data from Dokos et al. [7] well. However, the quality 
of fit was not equal for all shear directions. Specifically, in all trial runs the derived 
SEF did not fit to the shear data for the (nf) and (ns) direction, as shown in Figure 
1.8. As the magnitude of shear stress is much lower in those cases, it is likely that the 
fitness function neglects these shear modes. Previous studies have used coefficient 
of determination in the fitness function [18, 21], which will be a key consideration 
in our future work. 

To extend CHESRA to multiple experimental data sets, we included both shear 
[7] and biaxial data [9] in trial 2. Here, the biaxial data was fitted at the expense 
of the shear data. A possible cause for this could be a non-optimal individual and 
generation number. Since more data sets were given to the algorithm, optimization is 
challenging in only 50 generations. Therefore, increasing the number of generations 
and individuals may allow for better convergence and a more universal SEF to be 
found. 

Ensuring that CHESRA derives simple equations was important to this study. 
Therefore, we included the hyperparameter a to the fitness function to ensure penal- 
ization against long equations with nesting. While this is a valuable proof-of-concept, 
it has limitations. Specifically, it does not take computational time into consideration. 
Long but simple equations might be penalized more compared to shorter, but more 
complex or nested functions. The time it takes for an equation to be calculated may 
be a better measure of simplicity and will be explored in next steps. 
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Moreover, the current version of CHESRA can only generate SEF using quadratic 
and exponential functions; however, this could be widened to a larger function space. 
In addition, CHESRA may be extended to even more experimental data sets. Future 
work may also investigate to what extent the functions generated by CHESRA fulfill 
the mathematical requirements of SEF, such as convexity and strong ellipticity [11]. 
Lastly, it would be interesting to compare our method to a Multivariate Adaptive 
Regression Spline (MARS) model since studies by Jamei et al. [20] suggest that this 
can be more accurate than ESR. 


1.5 Conclusions 


In this work we present CHESRA, an ESR framework for automatically deriving 
strain-energy functions of the myocardium. Our algorithm allows for controlling 
function complexity and inputting multiple experimental data sources. 
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A 
Chapter 2 updates 
Electromechanical In Silico Testing Alters 
Predicted Drug-Induced Risk to Develop 
Torsade de Pointes 


Anna Busatto, Jonathan Krauß, Evianne Kruithof, Hermenegild Arevalo, and Ilse 
van Herck 


Abstract Torsade de Pointes (TdP) is a type of ventricular tachycardia that can occur 
as a side effect of several medications. The Comprehensive in vitro Proarrhythmia 
Assay (CiPA) is a novel testing paradigm that utilizes single cell electrophysiological 
simulations to predict TdP risk for drugs that could potentially be used clinically. 
However, the effects on mechanical performance and mechano-electrical feedback 
are neglected. Here, we demonstrate that including electromechanical simulations 
in CiPA testing can provide additional insights into the predicted drug-induced TdP 
risk. In this work, we analyzed six drugs, namely flecainide, ibutilide, metronida- 
zole, mexiletine, quinidine and ranolazine. We compared previously classified risks 
(low, intermediate, high) with our fully coupled electromechanical simulation results 
based upon the action potential, the electromechanical window, and the maximum 
active tension [1]. For ranolazine and metronidazole the predicted risk changed from 
low to intermediate and intermediate to high, respectively. For the latter, while elec- 
trophysiological markers indicated a low risk, the active tension decreased by 58% 
which can result in a loss of heart function. Therefore, adding mechanics to CiPA 
testing results in an altered prediction of drug-related TdP risk. 
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2.1 Introduction 


Torsade de Pointes (TdP) is a form of abnormal heart rhythm often preceded or 
caused by a prolonged QT-interval of the ECG [2]. Due to the high risk of unwanted 
and dangerous side-effects, promising new drugs with a potential TdP risk have been 
excluded from the research and drug-development pipeline. This, in combination 
with the long and complicated process of drug approval, led to a lack of new 
drugs on the market for cardiac disease. Therefore, the Comprehensive in vitro 
Proarrhythmia Assay (CiPA) is performed to better predict the TdP risk for drugs 
that could potentially be used in the clinic. However, some of the tests can result in 
inaccurate predictions of the actual outcomes and lead to restrictive results, limiting 
the number of drugs allowed on the market. 

There has been ongoing research on improving these models to increase the pre- 
dictability of drug-induced effects in the heart [3]. Llopis et al. proposed CiPA testing 
using a population of models approach to obtain more accurate risk predictions by 
simulating the effects of altered ion channel functions on the electrophysiological 
behavior of cardiac cells [1]. 

The electromechanical window (E Mw) is a biomarker for TdP risk introduced 
by Passini et al. for in silico drug testing [4]. Clinically, EM,, is defined as the 
difference between the duration of electrical and mechanical systole [4]. However, 
for the in silico test, they defined the EM, as the time difference between calcium 
transient duration at 90% repolarization (CaT Dog) and action potential duration at 
90% repolarization (APDo9) [4]. This assumption makes the simulation computa- 
tionally less expensive and complex, but this simplification reduces accuracy in the 
mechanical component. 

In previous work, the effects of electromechanical feedback are neglected and 
no insight into the mechanical function is obtained. This feedback alters the elec- 
trophysiological results, especially the calcium transient, and affects the biomarkers 
used to predict TdP risk. However, the effect of drugs on the mechanical performance 
can be investigated with electromechanical simulations. Therefore, in this work, a 
fully coupled electromechanical simulator SimCardEMS is used to perform CiPA 
tests [5]. Here, the cardiac mechanical function can be analyzed in terms of active 
tension, tissue deformations and EM,, where the latter can be altered due to elec- 
tromechanical feedback. We performed several tests to compare our predicted TdP 
risks to previous classifications such as the ones found in Llopis et al. [1]. 


2.2 Methods 


To assess TdP risk, we focused on selected biomarkers, namely maximum active 
tension Ta, max [kPa], APDoo [ms], CaT Doo [ms], and the EM, [ms]. Using these 
metrics, we analyzed six different drugs from CiPA: flecainide, ibutilide, metron- 
idazole, mexiletine, quinidine, and ranolazine with known cardiac effects. Each of 
the drugs has been classified as high, intermediate, or low risk for TdP in Llopis et 
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al. [1]. Here, the indicated risks are compared with the risks predicted by a fully 
coupled electromechanical model to investigate whether the mechanical biomarkers 
should also be taken into account for risk classification. 


2.2.1 Drug Implementation 


The cardiac effects of the drugs are tested for three different plasma concentrations 
(PC). For each of the concentrations, the ion channel conductances (g) are multiplied 
with scaling factors (SF) representing the effect of the tested drug. The SF values are 
calculated using the drug PC, the inhibitory concentration (C59) and Hill coefficient 
(h) value of the targeted ion channel: 


seer = et (2.1) 


SF = 
g ICso 


Table 2.1: Ion channel scaling factors (SF) for quinidine using three different plasma 
concentrations (PC). Scaling factors for seven ion channels were calculated with 
Equation 2.1. 


Scaling factor PCo.s PC\.0 PC29 
SFr 0.375 0.231 0.130 
S Fa 0.936 0.863 0.730 
SFyav 0.585 0.429 0.285 
S Feat 0.970 0.941 0.889 
SF xs 0.985 0.967 0.929 
SFxi 0.983 0.977 0.970 
SFio 0.731 0.524 0.309 


Equation 2.1 is derived from the standard Hill equation [1]. The PC, [C59 and h 
values are taken from Llopis et al.. The SF values are calculated for the considered 
ion channels Iķr, INa, Nat,» CaL, ks, [K1 and Žo at three different PCs which are half 
(PCo.s5), normal (PC;.9) and double (PC2.9) the effective free therapeutic plasma 
concentration [1]. As an example, specific SF values used to simulate the effect of 
quinidine are shown in Table 2.1. 


2.2.2 SimCardEMS Simulations 


The calculated SF values are finally used as input for the SimCardEMS solver, which 
is implemented in FEniCS [6], to simulate the effects of drugs on an endocardial tis- 


22 Electromechanical Testing in TdP 


sue block. We solved the electrophysiological O’ Hara-Rudy model which is strongly 
coupled with the mechanical Land model via intracellular calcium concentration and 
stretch activated channels [7, 8]. The active tension was scaled by a factor of five. 
The Holzapfel model was used to describe the tissue material properties, where we 
increased the stiffness by a factor of ten [9]. 

The simulations were performed on a 6x3x3 mm mesh representing an endocar- 
dial tissue slab with 1 mm and 0.5 mm spatial resolution for the mechanical and 
electrophysiological parts, respectively (Fig. 2.1). One corner of the mesh was fixed, 
and each of the connected planes were fixed in their respective normal directions, al- 
lowing the simulation of free contraction. The fibers were aligned in the longitudinal 
direction of the mesh and the tissue was activated by stimulating the entire mesh. 

To reach steady state, 40 beats of 1000 ms each were simulated and used as control. 
Thereafter, the drugs were added into the simulation by loading the calculated SFs. 
Each of these simulations was run either for five or ten beats (1000 ms each) 
depending on when steady state was reached; the goal was to investigate the effect 
of each tested drug on the electrophysiological and mechanical behavior of the cells. 
The biomarkers used for the analysis were extracted from the center of the mesh for 
the last beat of each simulation. 


2.2.3 Evaluation 


The TdP risk classifications were based upon maximum active tension Ta, max, 
APDo9, CaT Doo, and EM,, where 


EM, = CaT Do — AP Do0 (2.2) 


Fig. 2.1: (Left) Mechanical mesh 6x3x3 mm with 1 mm resolution. (Right) Elec- 
trophysiological mesh 6x3x3 mm with 0.5 mm resolution. Three planes are fixed in 
their normal direction to allow free contraction of the material. 
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The metrics used for evaluation are visualized in Fig. 2.2. TdP risk was classified as 
high or intermediate respectively in the following situations: 


e Ta max decreases by more than 50% (high) or between 20% and 50% (intermediate) 
compared to control. Clinically, a decrease in ejection fraction from 60% to 35% 
is considered dangerous [10]. We assumed that a 50% drop in Ta,max represented 
a similar level of danger for the patient since T} max affects the contraction; 

e APDo9 increases by more than 10% (high) or between 6% and 10% (intermediate) 
compared to control [11, 12]; 

e EM, decreases by more than 20% (high) or between 10% and 20% (intermediate) 
compared to control [4]. 


If none of these criteria were met, the drug was considered safe and classified as low 
risk. 
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Fig. 2.2: Transmembrane voltage, calcium concentration, and active tension for 


control. The metrics used to evaluate the results are indicated in the figure with 
dashed lines. 
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2.3 Results 


The biomarkers extracted for all simulations are shown in Table 2.2. For the control 
simulation we found a Ty max Of 55.0 kPa, APDoọo of 258 ms and EM of 134 ms. 
The associated risk classification thresholds for T4 max were 44.0 kPa (intermediate) 
and 27.5 kPa (high). For AP Dogo, the thresholds were 273 ms (intermediate) and 284 
ms (high), while for EM,, they were 121 ms (intermediate) and 107 ms (high). 


2.3.1 Action Potential Duration 


The resulting action potentials for control as well as the six examined drugs at PC1.0 
are shown in Fig. 2.3. Ibutilide, quinidine and flecainide have AP Doọ values of 597 
ms, 446 ms and 352 ms respectively which are above the threshold for high risk 
drugs. Ranolazine has an APDog of 283 ms at PC;.9 which fits in the range of 
intermediate risk drugs. APDog values for metronidazole and mexiletine are below 


Table 2.2: Overview of the extracted biomarkers, Ty max [kPa], APDoo [ms], CaT Doo 
[ms], EM,, [ms] and TdP risk classifications for all performed simulations. PCo.5, 
PCi. and PC3 o results are given in this order for each of the six tested drugs. 


Drug Tamaz AP Do CaTDọ%ə EM, in vitro in silico 
risk risk 

Control 55.0 258 392 134 - - 
61.0 316 404 88 

Flecainide 66.4 352 412 60 High High 
75.4 400 425 25 
50.8 531 541 9 

Ibutilide 46.1 597 560 3 High High 
43.8 644 641 -3 
28.8 249 393 145 

Metronidazole 22.9 252 396 144 Intermed High 
17.8 263 401 138 
54.8 253 390 136 

Mexiletine 54.3 251 388 137 Low Low 
53.2 250 386 136 
56.1 372 425 53 

Quinidine 57.3 421 446 25 High High 
61.9 470 469 -2 
68.7 272 391 119 

Ranolazine 79.5 283 391 107 Low Intermed 


96.1 303 392 89 
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the intermediate threshold and should therefore be classified as low risk if only the 
AP Dog is considered. 

All the chosen drugs were then examined for PCo.5 and PC2.9 as well. An example 
of the differences of action potential shape for quinidine based on the PC is shown in 
Fig. 2.4. A concentration dependent effect is observed for AP Dog as it is the smallest 
for PCo.5 and increases with an increase in PC. An overview of the APDoọ values 
for each of the examined drugs at PCo.5, PC1.0, and PC2.9 is given in Table 2.2. 


2.3.2 Electromechanical Window 


Looking at the EM,, values, high risk drugs were classified by 20% decrease com- 
pared to control; flecainide, ibutilide, and quinidine returned values below this thresh- 
old, indicating high risk, while metronidazole and mexiletine returned values above 
it. Finally, for ranolazine the decrease was slightly less than 20% which, in combi- 
nation with the remaining biomarkers, led us to classify it as an intermediate risk 
drug. For all the drugs, the TdP risk classification based on the EM,, was the same 
as APDop risk classifications. In four out of the six analyzed drugs, the E My, pro- 


control 
flecainide 
ibutilide 
metronidazole 
mexiletine 


quinidine 


ranolazine 


Vm in mV 


0 200 400 600 800 1000 
tin ms 


Fig. 2.3: Action potentials for control and six tested drugs at PC 9. Based on these 
results, ibutilide, quinidine and flecainide are classified as high risk. Ranolazine is 
classified as intermediate risk for this drug concentration. 
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gressively decreased with higher drug concentration compared to control. However, 
addition of metronidazole at PCo.5 increased the EM,, initially while for higher 
concentrations of metronidazole (PC2,9) the EM,, decreased again. The addition 
of mexiletine changed the EM,, by only 3 ms while both the APDoo and CaT Doo 
remained similar to control. All results are shown in Table 2.2. 


2.3.3 Maximum Active Tension 


Applying metronidazole, the maximum active tension decreased from 55.0 kPa to 
28.8 kPa, 22.9 kPa, and 17.8 kPa for PCo.5, PCi.9 and PC2.9, respectively. This 
represents a decrease in Ty max Of 48%, 58%, and 68%. For ranolazine, the Ta max 
increased by 45% at PC; 9 compared to control. In the remaining drugs, the maximum 
active tension remained similar or slightly increased compared to control. 


control 


=== quinidine at PCo.5 
quinidine at PC, 0 


quinidine at PC2 9 


Vm in mV 


0 200 400 600 800 1000 
t in ms 


Fig. 2.4: Action potential for control and quinidine at PCo.5, PC1.0, and PC2 0. 
An increase in PC highlights the drug concentration effect on the action potential 
duration. 
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2.3.4 In Silico and in Vitro TdP Risk Classifications 


We classified all six drugs using the evaluation process described in section 2.2.3. 
Flecainide, ibutilide, and quinidine were determined to be high risk both in the 
previous classification and in our electromechanical model (based on APDop and 
EM,, values). Based on the maximum active tension, the risk of metronidazole 
was changed from intermediate (CiPA indication) to high risk. On the other hand, 
following the APD and EM,, values it would have been classified as low risk. 
Mexiletine was classified as low in the CiPA as well as in our study. Lastly, ranolazine 
was previously predicted as low risk but was classified as intermediate risk using 
our model based on both the APDoo and E Mw values. Even though the increase in 
Ta,max for ranolazine was the largest of the analyzed drugs, potentially having an 
effect on contraction and tissue stress, increases in Ta max Were not considered in TdP 
risk classification. 


2.4 Conclusions 


In this study we analyzed six drugs, flecainide, ibutilide, metronidazole, mexiletine, 
quinidine, and ranolazine. Three of these drugs were classified as high risk, one as 
intermediate risk, and two as low risk according to CiPA. Using our fully coupled 
electromechanics model, four drugs were classified as high risk, one as intermediate 
risk and one as low risk. An overview of this distribution is given in Table 2.3. After 
performing our analysis, we concluded that some drugs may belong to different 
risk categories depending on the individual parameters considered. For example, 
ranolazine was previously classified as low risk; in our simulation, both the AP Doo 
and E'M,, biomarker values for PC; fell in the range of intermediate risk, while for 
PC?,9 it was predicted high risk. In general, for all six drugs, both the AP Doo and the 
EM,, biomarkers indicated the same risk categories. However, the risk classification 
based on Ty max differed from the other biomarkers. For metronidazole, a low TdP 
risk was observed based on the electrophysiological biomarkers. Despite that, there 
was a decrease of 58% in Ta,max Which can lead to severe problems in heart function. 


Table 2.3: Predicted versus in vitro risk classifications for the analyzed drugs. Out 
of the six tested drugs, two drugs were classified differently based on the electrome- 
chanical results. 


Pred. high Pred. intermediate Pred. low 
CiPA high 3 0 0 
CiPA intermediate 1 0 0 


CiPA low 0 1 1 
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We conclude that the addition of a mechanical biomarker, the maximum active 
tension, is valuable in the prediction of cardiac risk stratification. Additionally, using 
the electrophysiological biomarkers, SimCardEMS succeeded in predicting the same 
TdP risk category as previously determined for four of the six drugs. 
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Osteoarthritis in Human Articular 
Chondrocytes 
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Maleckar 


Abstract Osteoarthritis (OA), a progressive degenerative disease of cartilage in 
joints, is the most common cause of chronic disability in older adults. While OA 
is mostly considered an age-related pathology, women have a 1.5-fold higher risk 
of developing OA relative to men and experience more severe symptoms. Yet, they 
remain underrepresented in musculoskeletal research and clinical trials. Responsible 
for cartilage formation, articular chondrocytes experience physiological changes in 
OA, but the functional implications of such alterations remain largely unexplored due 
to difficulties in acquiring the data experimentally. Through reparameterization, we 
expand a mathematical chondrocyte model to investigate sex-specific OA pathogen- 
esis. We performed sensitivity analysis to address the impact of ion channel activity 
in healthy and OA chondrocyte populations. Simulations show that in healthy female 
chondrocytes, the resting membrane potential is more depolarized than in healthy 
male chondrocytes, suggesting potential sex-specific emergent physiological differ- 
ences in articular chondrocytes. In both sexes, the resting membrane potential of 
healthy chondrocytes is most sensitive to ICa-ATP, INa-b, [Nak and Iķg-p, but in 
OA it depolarizes and becomes sensitive to IkpR, Inax and Ig-p. Developed and 
evaluated against experimental data, our articular chondrocyte OA electrophysiolog- 
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ical model can be used to further study OA pathology and sex-specific pathological 
OA changes. 


3.1 Introduction 


Mammalian chondrocytes are found within intervertebral discs and articular carti- 
lage, where they maintain the extracellular matrix (ECM) and produce the cartilage 
matrix [1]. The ECM mainly consists of collagen fibers that transmit force and dis- 
sipate energy, proteoglycans that withstand compressional forces, and synovial fluid 
to reduce friction between bones. In normal function, healthy chondrocytes respond 
to outside stimuli and damaged articular cartilage, where they increase production 
of specific ECM complexes [2]. Functional, mature, healthy articular chondrocytes 
are necessary for the maintenance of cartilage and therefore, physiological joint 
motion [3]. Chondrocytes function to maintain a homeostatic environment in the 
articular cartilage. Pathophysiologic conditions can lead to the development of os- 
teoarthritis (OA); however, the mechanism by which this occurs is poorly understood. 
Osteoarthritis is the most common form of arthritis and is a disease characterized 
by the degeneration of cartilage and underlying bone [4]. Women are 1.5x more 
likely to develop OA, and experience more severe symptoms as compared to men 
[5, 6, 7]. In OA, chondrocytes secrete increased levels of inflammatory cytokines, 
actively produce proteoglycans and collagen type II to recover the degeneration of 
the ECM, and become hypertrophic [8]. The present study on human, sex-specific 
OA is modeled after the specific cell physiology of the human articular chondrocyte. 

Recent research suggests that maintaining a stable resting membrane potential 
is essential for chondrocytes in articular cartilage to be able to withstand pressure 
and force changes [9]. Despite chondrocytes being non-excitable cells, they undergo 
robust ion-channel mediated changes in response to OA. This is further supported 
by experiments where ion-channels contributing to the resting membrane potential 
were blocked, resulting in a decrease in the production of matrix mRNAs, proteins 
and glycosaminoglycans [10, 11]. Additionally, channel blockers inhibit chondro- 
cyte proliferation and increase apoptosis [12, 13]. During OA, chondrocytes have 
reduced capability to respond to changes in the extracellular environment, specifi- 
cally with respect to ion transport; this is recognized by i.e. a deficiency in volume 
regulation [14]. Chondrocyte diameters vary between 7 and 30 pm, making technical 
electrophysiological experimental studies difficult [15]. Therefore, there is limited 
experimental data investigating the electrophysiological differences of chondrocytes 
in normal physiology and pathophysiology. 

Therefore, to gain a better understanding of the human chondrocyte channelome in 
the context of OA, we have developed an OA chondrocyte model by further expanding 
a previously established mathematical model of the human articular chondrocyte 
[3, 16]. Mathematical models can play a key role in the dearth of physical data and 
in hypothesis generation, among other utilities. In the absence of constraining data, 
models can be used to generate valuable hypotheses about the underlying processes 
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at work, which can then be tested in carefully designed experiments and/or against 
data from other sources. This initial model of the human articular chondrocyte 
features various ion channels describing the resting membrane potential and ion 
movement across the cell membrane. Previous iterations of the model focused on 
potassium fluxes across the cell membrane, and included functions specifically for the 
integration of the sodium/potassium pump, an active transport pump that establishes 
an electrochemical gradient across the cell membrane. In this study, the model 
is expanded to include sex-specific OA by implementing the differences in ion 
channel conductances observed between males and females, as well as alterations 
of ion channel expression in OA disease. Populations of models using a log-normal 
distribution were also implemented to accurately represent biological variability, 
and subsequent parameter sensitivity analysis was performed to identify the currents 
that influence chondrocyte resting membrane potential to the greatest degree in both 
health and disease. 

In summary, we have developed a mathematical OA chondrocyte model to gain 
a better understanding of OA pathogenesis through exploratory simulations, driving 
forward hypothesis development in the context of relative scarcity of human chon- 
drocyte electrophysiological data across the channelome. As OA is more prevalent 
in females than males, we also implemented sex differences in our updated model. 
Our main goals in the present study were to: (i) update and scale the baseline male 
model (1) to generate a (2) disease OA-male model, (3) a control female model, 
and (4) a disease OA-female model, (ii) utilize a population of models approach 
to determine which ion channels were the main contributors to resting membrane 
potential in human articular chondrocytes during normal function and in OA. 


3.2 Methods 


To study the effects of OA on the resting membrane potential, we employed the 
chondrocyte model and simulation protocols developed by Maleckar et al. [3] and 
further improved by Fischer-Holzhausen et al. [17]. The model components are 
shown in Figure 3.1. Iyag has been scaled to be within 0.1-0.6 of the original value, 
corresponding roughly to 0.9-5.35 pA/pF, as prior work has shown this range is 
more likely to present a physiologically-relevant model regime [17]. 

Following the Hodgkin-Huxley formalism, the cell membrane is modeled as a 
capacitor coupled in parallel with ion channels represented by resistors. The time- 
evolution of the transmembrane potential is given by the summation of all electro- 
genic transport processes. 


dV 
Cmr = —(IkDR + [Nak + lNaCa + ICa-ATP +1K-ATP (3.1) 


+ Tk2pore t+ nab tl xp +Icib +Ikca) , 
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Fig. 3.1: Schematic illustration of the model adapted from Maleckar et al. [17], 
showing the principal ion-selective channels, exchangers, and pumps expressed in 
human chondrocytes. 


where Cm is the chondrocyte capacitance (6.3 pF). A set of three transport equations 
gives the time-derivatives of the intracellular ionic concentrations of Nat, K*, and 
C at 


d[Na* ]; _ _INab+3INak +31Naca-INan + Icıb 


= ; 3.2 
dt voli X F 02) 
d[K*]; Txp -2Inak +IKDR+ IK2pore + IKCa+IK-ATP 
== : (3.3) 
dt voli X F 

d |Ca”]. Ica- Fie, dO. 
Lea]: _ _tca-are —2Inac 0.045 x = , (3.4) 

dt 2vol, X F dt 


where vol; is the internal volume of the chondrocyte (0.005884 mL), F is the Faraday 
constant (96,485 C/mol), and Oc is the fraction of intracellular calmodulin bound 
to Ca”. The coefficients in front of the exchangers come from their respective ionic 
exchange ratios. For example, the NaK exchangers exchange three Na* ions for two 
K* ions. 
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3.2.1 Modeling Ionic Changes Induced by Osteoarthritis in Male and 
Female Chondrocytes 


To gain mechanistic insight into the pathogenesis of OA, we implemented changes in 
ion channel expressions as measured by Karlsson et al. [18]. We further expanded our 
investigation by integrating sex-specific differences in ion channel activity into our 
model. Due to the lack of experimental data on female chondrocytes, we compiled 
a list of sex-specific differences in ion channel expression based on mRNA and 
ion channel subunit expression data collected from epicardial cardiomyocytes. In 
our model, we modeled the effect of changes in ion channel expression on their 
conductances (G), as displayed in Table 7.2. 


Table 3.1: Sex-specific subcellular ion channel conductances relative to the male 
baseline model 


Parameter Male Control Male OA Female Control Female OA 

GNaK 1 2.20 0.95 2.08 

GCa-ATP 1 2 1 2 

GK-2Pore 1 0.20 1 0.20 

GKDR 1 8.30 0.80 6.64 i 
GK-ATP“ 1 0.34 1 0.34 

GNaCa 1 1 0.97 1 

GK-b 1 1 0.51 1 


Implemented as change to the Q10 parameter which defines the effect of temperature on ion channel 
kinetics [19, 20]. 


3.2.2 Generation of a Population of Models 


To capture biological variability and investigate parameter sensitivity, we used the 
published chondrocyte model [3, 17] as a baseline to build our population of 1,000 
models by randomly modifying parameters corresponding to conductances of ion 
current and maximal rates of ion transports. A relatively large sample size of 1,000, 
typically used in similar population-based studies [21, 22], was chosen to allow 
even small effects, often indistinguishable in small sample sizes, to reach statistical 
significance. Prior studies have shown that biological parameters, when randomly 
selected, are distributed according to log-normal frequency curves instead of normal 
ones [23]. For each channel in the model, we varied its conductance by sampling 
from a log-normal distribution centered around their respective baseline value with 
the standard deviation set to 0.15. 
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3.2.3 Model Parameter Sensitivity Analysis 


We performed parameter sensitivity analysis on the generated model population to 
assess the sensitivity of the model output (e.g., the resting membrane potential) in 
response to channel conductance variations. The degree to which perturbing a set 
parameter can influence the results are quantified as a set of regression coefficients, 
which are obtained by performing multivariable partial least squares regression on 
the dataset. Partial least squares regression was chosen over standard multivariable 
regression as many independent parameters (e.g., ion channel conductances) were 
used to predict a smaller set of dependent variables (e.g., the resting membrane 
voltage) [24]. This methodology was first used in cardiac electrophysiology by 
Sobie et al. [25] and has been widely used in other fields of biology [21, 26]. 


3.3 Results 


3.3.1 Modeling the Impact of Sex-Specific OA on Resting Membrane 
Potential 


In Figure 3.2, we implemented the experimentally-observed changes in ion channel 
expression in control and OA chondrocytes, accounting for sex differences as de- 
scribed in Table 7.2. Figure 3.2 shows that the healthy female chondrocyte maintains 
a more depolarized resting membrane potential (-58.21 + 5.75 mV) as compared to 
the healthy male chondrocyte (-69.11 + 4.71 mV). Electrophysiological differences 
instigated by OA result in depolarization in both male (from -69.11 + 4.71 mV to 
-53.87 + 3.21 mV) and female (from -58.21 + 5.75 mV to -49.03 + 1.34 mV) chon- 
drocytes. This overall depolarizing effect of OA remodeling on the resting membrane 
potential has previously been reported in experiments [27]. 
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Fig. 3.2: Resting membrane potential in control and OA chondrocytes for (A) male 
and (B) female models. Changes in ion channel expression in OA induce resting 
membrane depolarization. 
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3.3.2 Population of Models 


Shown in Figure 3.3, and as outlined previously, we further expanded our analysis to 
investigate population variability by generating 1,000 models for each condition with 
channel conductances varied according to a log-normal distribution. For all 1,000 
parameter sets, the steady-state solutions were reached within the simulation period 
of 50,000 seconds. Figure 3.3 shows that in male and female OA chondrocytes, 
intracellular sodium quickly becomes depleted and remains very small, close to 
zero, for the rest of the simulation. We calculated the mean resting membrane 
potentials for all conditions and compared them against the reported experimental 
values as measured in experiments, shown in Figure 3.4. In all cases, the membrane 
potentials of OA chondrocytes are relatively more depolarized compared to their 
healthy counterparts. Comparing male and female chondrocytes, although the control 
resting membrane potential values slightly differ (females are more depolarized 
compared to males), OA male and female chondrocytes depolarize to a similar value 
(approximately -51 mV). 


3.3.3 Parameter Sensitivity Analysis 


Sensitivity analyses were performed on the chondrocyte populations to reveal how 
ionic modulations induced by OA might affect the sensitivities of underlying elec- 
trophysiological properties in the models. Shown in Figure 3.5, in the healthy male 
chondrocyte population, Ij ax, Ica-aTP> INa-p, and Ix _p are predicted to have the 
most critical impact on the resting membrane potential. In this case, Iya-» is the 
strongest depolarizing current and Jcqg_ arp is the strongest hyperpolarizing current. 
In male chondrocytes expressing OA, only Iyax and Ig-p remain the most influen- 
tial currents. Compared to the healthy case, the impact of Icg—arp and Iya-p are 
severely reduced in OA, while Iyax gains a 1.6-fold depolarizing influence. 

The healthy female chondrocyte population shares some similarities with their 
male counterparts, as the resting membrane potential is still the most sensitive to 
INaK, ICa-ATP, INa-b, and Ig-p. However, Ig-arp now plays a slightly more 
dominant role as a hyperpolarizing current. Compared to male chondrocytes, Inax 
gains a more significant role in establishing the resting membrane potential than in 
males (2.5-fold increase in importance). Iyag is the strongest depolarizing current 
and [x~» is the strongest hyperpolarizing current. Surprisingly, in female chondro- 
cytes expressing OA, the majority of the currents no longer have as much of an 
impact on the resting membrane potential. Nonetheless, Iyag and Ig-p each still 
retain a minor influence. 
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Fig. 3.3: Membrane voltage and internal ion concentrations for (A) male and (B) 
female population. 


3.3.4 Inhibiting Iyag Restores Normal Resting Membrane Potential in 
OA Chondrocytes 


Per sensitivity analysis results, Iygx has the largest impact on depolarizing the 
resting membrane potential (i.e., having positive regression coefficients) in both 
male and female chondrocytes in OA. To restore the membrane potential back to 
near control values, we applied jax current block treatment then simulated the 
membrane potential in male and female OA populations of 1,000 models each. Figure 
3.6 shows that applying 45% Inax block to male and 55% Iyax block to female 
OA chondrocytes restored the mean resting membrane potential of the population 
back to the control level. Notably, previous clinical studies have revealed that Inax 
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Fig. 3.4: Mean resting membrane potential in control and OA chondrocytes for male 
and female populations (1,000 models each) compared against experimental data 
[27]. To account for differences in data collection methodology, experimental resting 
membrane potential values were scaled to match the mean of the male population. 
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Fig. 3.5: Parameter sensitivity analysis on (A) male and (B) female population of 
sex-specific chondrocyte models. The regression coefficients measure the relative 
impact of a particular ionic current on the resting membrane potential. A positive 
regression coefficient indicates the associated current is a depolarizing current, while 
a negative regression coefficient indicates a hyperpolarizing current. 
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blockers Ouabain and Digoxin protected against OA and relieved OA-associated pain 
[28]. 
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Fig. 3.6: Resting membrane potential in control, OA, and OA chondrocytes with 
partial block of Iyax for male and female model population (1,000 models each). 
Blocking Iyax (45% for male and 55% for female population) helped restore the 
resting membrane potential back to normal levels in OA chondrocytes. 


3.4 Discussion 


OA, a debilitating disorder affecting 1 in 7 individuals in the United States in their 
lifetime [29], has no current curative treatment. The prevalence, as well as disease 
severity, is more pronounced in female patients [30]. It is known that chondrocytes 
play an essential role in maintaining healthy cartilage, and that their function is 
impaired in disease [31]. As experimental investigation of human chondrocytes is 
challenging due to the lack of human control samples, we have further expanded 
a previously-developed mathematical model of a chondrocyte channelome by in- 
troducing experimental changes observed in ion channel electrophysiology, as well 
as temperature and capacitance [3]. The updated model is user-friendly and freely 
accessible (github. com/k-ngo/Chondrocyte). 

Simulations revealed that the chondrocyte membrane becomes depolarized in 
OA, which agrees with experimental data [32]. This suggests that the changes in 
ion channel expression together with temperature and capacitance are sufficient to 
capture the main pathological alterations observed in OA chondrocytes related to 
their membrane potential. In the future, it remains important to further validate the 
model experimentally, as well as to investigate the primary cause of the alterations 
in ion channel expression observed in OA chondrocytes. These are likely to involve 
inflammatory processes associated with OA [33]. 
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Comparisons between male and female chondrocyte models reveal intriguing 
phenotypic differences. First, the resting membrane potential is more depolarized 
in females compared to males in controls. The main differences between male and 
female healthy chondrocyte models include reduced Gxpr (by 25%), Gx_p (by 
95.2%), G NaCa (by 2.8%) and G yag (by 5.7%). The investigation of the main fac- 
tors important for maintaining membrane potential visualized in Figure 3.5 revealed 
that, in the present model, calcium ATPase and background sodium and potassium 
currents contribute the most for maintaining healthy resting membrane potential in 
males, with an addition of Jyax in females. It could be that the relative depolariza- 
tion of the female chondrocyte membrane driven by the aforementioned changes in 
the model predisposes female patients to OA, although additional experimental data 
is needed to further investigate this hypothesis. 

OA chondrocyte model simulations show that resting membrane potential is de- 
polarized in disease, following experimental data [32]. This effect is seen in both 
male and female OA models. Our results also suggest that chondrocyte resting mem- 
brane potential in female patients could be more depolarized compared to that in 
males, which could lead to higher depolarization in OA (mean resting membrane 
potential is -49.03 + 1.34 mV in the female and -53.87 + 3.21 mV in the male OA 
model). The resting membrane potential in OA was around 17.5% more depolarized 
for female and 24.5% for male model compared to respective healthy resting mem- 
brane potentials. In our OA model, yax current plays a prominent role in regulating 
membrane voltage: in both males and females with OA, increased Iyag becomes the 
most important contributing factor to the overall depolarization of the resting mem- 
brane potential, as revealed by sensitivity analysis. Inhibiting Iyag might therefore 
provide a useful approach for reversing resting membrane alterations in OA and, 
more broadly, Jax could be a useful target in modulating disease pathogenesis in 
OA. 

Tnak inhibitors are well characterized for e.g. cardiac disease treatment, and new 
selective inhibitor drugs are currently in the pipeline [28]. To test whether Iyax 
block could be useful in overcoming OA-associated changes in resting membrane 
potential, we simulated Zyag block. With varying degrees of Iyax inhibition, the 
resting membrane potential in disease can be restored to healthy, pre-OA levels in both 
males and females. Therefore, re-purposing [jax blockers feasibly provide a useful 
new therapeutic approach for OA treatment. [ax blockers have, in fact, already 
been applied for OA treatment: clinical studies have reported positive outcomes in 
patients, revealing that Iyag blockers, Ouabain and Digoxin, induced a reduction in 
pain and elicited an apparent OA-protective effect [28], supporting our simulation 
results implying J)ax blockers as a potential treatment for OA. 

Given this putative critical role of Na+ concentration in healthy and OA chon- 
drocytes, novel future experiments should consider thoroughly probing Na+ balance 
in human chondrocytes. For instance, Figure 3.3 shows that in male and female 
OA chondrocytes, intracellular sodium quickly becomes depleted and remains very 
small, close to zero, for the rest of the simulation, suggesting that alternate sodium 
regulation not currently accounted for in base models may apply in OA. Model 
development should additionally move towards incorporating new channel species 
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to regulate Na+ concentrations in control and in OA, including epithelial sodium 
channels (ENaC), which are not currently included in the model [31, 34]. Other work 
has also focused on voltage-gated sodium channels in chondrocytes, which may also 
be considered. 

Naturally, the lack of available biological data from new electrophysiological and 
other experimental studies presents a challenge for the validation of both male and 
female models, both in control and in OA. Further experimental studies will be 
instrumental in determining the ranges of parameters for the model, as well as for 
simulated chondrocyte behavior evaluation. In the future, further experimental data 
will be helpful in further developing and expanding the current chondrocyte model. 
While the current study, synthetic in nature, limits its direct translation to clinical 
application, the present work retains utility in terms of physiological exploration and 
hypothesis generation. 

In summary, further incorporating additional species key in the chondrocyte 
channelome into the model and incorporating new data from emerging experiments 
will permit even greater insights into the broad and nuanced role of membrane 
dynamics in chondrocyte function and pathophysiology. 


3.5 Conclusions 


We have further expanded a previously published chondrocyte mathematical model 
enabling simulations of the osteoarthritic chondrocyte channelome. Implemented 
changes were based on published data, and model simulation results show that the 
resting membrane potential in chondrocytes becomes depolarized in OA in both 
male and female chondrocytes, which agrees with the experimental data acquired 
from males. Our sex-specific chondrocyte model reveals differences in resting mem- 
brane potential between males and females that could potentially contribute to the 
higher prevalence of OA in female patients. Finally, sensitivity analyses revealed the 
main currents responsible for maintaining resting membrane potential in chondro- 
cytes and potential novel therapeutic targets for OA treatment. Our OA chondrocyte 
electrophysiological model provides an accessible tool for subsequent studies of OA 
pathogenesis and drug targeting. 
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A 
Chapter 4 updaies | 
Recapitulating Functional Heterogeneity in 
Electrophysiologically Active Tissues 


Meye Bloothooft, Joseph G Shuttleworth, Gabriel Neiman, Ishan Goswami, and 
Andrew G Edwards 


Abstract Inter-cellular heterogeneity is central to the dynamic range and robustness 
of function in many tissues, particularly electrically excitable tissues. In pancreatic 
islet 6-cells, inter-cellular heterogeneity underlies the range of insulin response to 
glucose. In human-induced-pluripotent stem cell-derived cardiomyocytes (hiPSC- 
CMs), inter-cellular heterogeneity presents a key challenge for drug screening appli- 
cations. In this study, we assess the ability to reconstruct inter-cellular heterogeneity 
in silico by applying a “population of models” (PoMs) framework, i.e. collections 
of computational cells created via Monte Carlo variation of model parameters. We 
define parameter variation based on experimentally observed heterogeneity in prop- 
erties such as ion current conductances and enzymatic affinities. We then assess the 
accuracy of those reconstructions, based on the degree to which variation in POM 
outputs (e.g. action potential duration) matches experimentally observed variation. 
We report that this “ground-up” approach underestimates functional heterogeneity 
in the hiPSC-CM population, but overestimates it in adult human cardiomyocytes. 
In contrast, the 6-cell PoM captures three distinct and physiologically relevant sub- 
classes of 6-cell function. In the future, we expect PoM approaches like these will 
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permit incorporation of realistic cellular heterogeneity in detailed models of intact 
tissues, and thereby aid development of sophisticated tissue-engineered platforms 
for therapeutics. 


4.1 Introduction 


Functional heterogeneity among cells within a tissue is a characteristic and relevant 
feature of living systems. Even within a genetically homogeneous population, in- 
dividual cells often show a considerable amount of phenotypic heterogeneity [1], 
and this feature has implications in the plasticity and robustness of the organisms. 
Plasticity allows adaptive cellular responses to perturbations in the environment 
while maintaining a robustness in functional behaviour of the tissue even after a 
perturbation or insult. For example, when a region of the heart becomes infarcted 
(dies), non-infarcted healthy areas will adapt and increase their force of contraction 
to maintain blood flow to the body. This ability to adapt is a form of robustness of 
organ function, but must be achieved in a manner that maintains the heart’s elec- 
trophysiologic stability. Similarly, in obesity induced diabetes, hyperlipidemia and 
hyperglycemia exert excessive insulin demand on pancreatic islet B-cells. In some 
B-cells this prolonged exertion leads to exhaustion of insulin reserve and cell death. 
Thus, increased demand for insulin must be achieved via an adaptive response in 
the remaining viable B-cells, which are phenotypically distinct. In both of these 
examples, cellular heterogeneity is central to the plasticity and robustness that al- 
lows both organs to maintain function after pathophysiological challenge. Functional 
heterogeneity may also broaden the range of responses to drugs [2]. For example, 
glibenclamide has differential effects on pancreatic islets owing to functional het- 
erogeneity of their B-cell populations [3]. Thus, the study of cellular heterogeneity 
for both basic understanding of living systems and discovery of therapeutics has 
motivated the development of experimental techniques and the construction of com- 
putational models to simulate the dynamics and effects of cellular heterogeneity 
[4]. 

A range of single-cell experimental techniques are available to investigate this 
heterogeneity. These range from transcriptome assessment via single cell RNA se- 
quencing, to classical cell electrophysiology via patch clamp. While these approaches 
are generally both time and labor intensive, data collected by those methods can be 
stored and interrogated via efficient modelling frameworks, such as the population 
of models approach (PoM) [5, 6]. These approaches simulate the dynamics and 
functional implications of heterogeneity in cell populations, including responses to 
pathologic insults and drug challenges. 

In the PoM approach, heterogeneity is implemented by varying parameters (such 
as the ionic conductance of a single or multiple ion channel/s) of a baseline model, 
which itself captures the average characteristics of a particular cell type. Each unique 
configuration of parameters creates a new instance of the cell model, which is anal- 
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ogous to anew “member” cell of the population. In general, this procedure has been 
performed by first specifying equal distributions over which all parameters of inter- 
est (e.g. all ionic conductances) can vary, and then selecting only members of that 
population that exhibit behaviour within the experimentally observable range [5, 6]. 
For example, a cardiac myocyte population could be constructed by first allowing 
all ionic conductance parameters to be sampled from Gaussian distributions with an 
arbitrary standard deviation of 30% of each mean parameter value in the baseline 
model - this sampling strategy involves homogeneous variability across parameters. 
Some members of this PoM will exhibit parameter combinations that are biolog- 
ically unrealistic, and thus possibly exhibit behaviour outside the experimentally 
observable range. To eliminate those models, it is common to reject models that 
produce outputs or features, such as a cardiomyocyte action potential duration, that 
are outside the range of experimental observations for the respective features. Con- 
tinuing the example, a cardiomyocyte population member will be deemed non-viable 
if the action potential duration is far longer or shorter than observed experimentally. 
Similarly, a non-firing B-cell model will be rejected from the islet PoM. Because, 
in these examples, the PoMs are generated from arbitrary and homogeneous pa- 
rameter distributions, they will exhibit appropriate behavioural characteristics, but 
may not involve truly biologically representative parameter combinations. Thus, the 
behaviours may not result from truly representative underlying dynamics. 

In this work, we apply a different approach, with the intention of introducing more 
biologically realistic variability in PoMs of three different cell types: human adult 
cardiomyocytes, hiPSC-CMs and primary islet 6-cells. By introducing variability to 
the input parameters that is grounded in experimentally observed variability for each 
parameter, we hope to span the biologically reachable configurations and limit the 
number of unrealistic configurations that are generated. That is, we use published 
experimental variation (standard deviations) for ionic conductance parameters and 
enzymatic activity/affinity parameters in order to apply parameter-specific hetero- 
geneity rather than variation that is homogeneous across the varied parameters. 
Based on these constraining conditions, we demonstrate the feasibility of PoMs to 
generate a heterogeneous population representing different functional phenotypes. 


4.2 Methods 


An illustrative example of our approach in generating PoMs is shown in Figure 
4.1. To recapitulate the functional heterogeneity in the tissue PoMs, we introduced 
parameter-specific variability into our models. To do this, we identified a set of pa- 
rameters which we expect to account for much of the inherent biological variability 
influencing electrophysiologic function in these cell types. These parameters are a 
selection of maximal conductance (or flux) parameters for both the cardiomyocyte 
models (hiPSC-CM and adult) and the 6-cell model. Additionally for the 6-cell, a 
selection of kinetic parameters for enzyme reactions in the glycolysis cycle were 
varied to introduce metabolic heterogeneity. We chose to endow each of these pa- 
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Fig. 4.1: Illustrative example of generating population of model (PoM) of human 
cardiomyocytes and islet B-cells. Experimental data obtained from literature was 
used to generate log-normal distributions of electrophysiological and/or metabolic 
parameters used as inputs to computational model. Sets of features are obtained 
from the simulations that are then used to generate viable candidates for the PoMs 
representative of cardiomyocyte or 6-cells via constraining conditions. 
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rameters, p, with a log-normal distribution, such that log p ~ pN(0,o7) where p 
is the original, unmodified parameter from the baseline model and ø is a standard 


deviation derived from the literature. We term the random variable P the scaling 


factor. To assign values for the various o, we surveyed the literature for measured 
variability of each parameter of interest. To mitigate the impact of systematic differ- 
ences among the means of different datasets, and to keep our parameter distributions 
centred on the original parameter values, we computed the coefficient of variation 
for each parameter (cp): 
ae 
Hiit 
where ujit and Ciit are the mean and standard deviation, respectively, of the observa- 
tions found in the literature. 

The published mathematical models chosen as the baselines for the adult-CM and 
the hiPSC-CM were developed by [7] and [8], respectively. For the 6-cell PoM, the 
electrophysiological model developed in [9] and glycolysis model developed in [10] 
were used to generate the heterogeneous population of f-cells. In the following sub- 


sections, we will detail the generation of parametric distributions and the constraining 
conditions via feature extraction for each PoM. 


Tp: 


(4.1) 


4.2.1 Cardiomyocytes 


The O’ Hara-Rudy model (ORd, [7]) is a widely-used mathematical description of the 
average human left ventricular cardiomyocyte. We specifically used the epicardial 
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version for our adult human baseline model. Likewise, the baseline model for hiPSC- 
CMs was that published in [8]. Both models use the following fundamental equation 
for calculation of time-varying membrane potential: 


dV 
nY 4.2 
Cm 2, (4.2) 


where C,,, is the capacitance of the cell membrane and each 7; is a current describing 
the transmembrane current carried by each species of ion channel present in the 
model. For each ion channel, there is a maximal conductance parameter g;, or analo- 
gous permeability parameter, such that J; = g; Î for some /;. /; is in turn described by 
a system of ordinary differential equations that specify the voltage, ion concentration 
and time-dependencies of ion channel function. For all currents, we only vary the 
maximal conductance (g;), and not the parameters of the remaining equations con- 
tributing to Î;. The formulation for sarcoplasmic reticulum (SR) calcium reuptake is 
different to those of the ion channels, but an analogous term that reflects the maxi- 
mal rate of SR calcium reuptake Jserca. As for the ion channel conductances, this 
parameter was varied to introduce heterogeneity in this key element of intracellular 
calcium handling. 


Table 4.1: Scaling parameters for the cardiomyocyte PoM 
Sources to construct coefficient of variation (op) listed by the reference numbers 


Variable Mean Unit op 
human adult CM 

&Kl 3.60 pA/pF 0.78! 

2kr 0.31 pA/pF 0.13! 

8&Ks 0.18 pA/pF 0.657 

ko 9.30 pA/pF 0.99! 
2CaL 10.20 pA/pF 0.13! 

2Na 16.10 pA/pF 0.837 
&NaL 0.34 pA/pF 0.494 
2NaK 1.90 jM/min/mg protein 0.485 
2NaCa 0.05 nmol Ca2+/mg protein/sec 0.3667 
JSERCA 7.10 nmol ATP/mg protein/min 0.395 
hiPSC-CM 

2Na 166.00 pA/pF 0.2681! 
8Ca 12.20 pA/pF 0.41 811-15 
kr 2.02 pA/pF 0.478:11,16-18 
&Ks 1.30 pA/pF 0.58811.19 
gr 2.50 pA/pF 0.32811 
2k1 1.00 pA/pF 2.2081 1,17,20 
gio 4.69 pA/pF 0.488-14-21 
8&NCX 2.10 pA/pF 0.3322 
2NaL 0.70 pA/pF 0.86!° 


H11], [12], °[13], 4[14], 5015], [16] 7[17], 3118], 9119}, '°(20), [21], 17[22], 7123], “24, 
15[25], 16/26], 17027], 18[28], 19/99], 20130], 21/731], 22132] 
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Fig. 4.2: A. The resultant log-normal distributions for the gx; scaling factor, along 
with histograms resulting from 100 samples, are shown for the representative para- 
metric variation in adult CM and the hiPSC-CM inputs for the models. B. The scaling 
factors that were sampled in order to create the PoM are shown for two cardiomyocyte 
models. Parameters from the models that were excluded from the POM are shown in 
grey, whereas the parameters which were calibrated into the PoMs are highlighted 
in orange. 


Table 4.1 provides the values for each varied parameter for both the adult human 
myocyte model and the hiPSC-CM model. The variability presented in this table 
may not completely capture the variability in each maximal conductance parameter. 
Nevertheless, these values provide a reasonable approximation, and yield plausible 
distributions for each of the relevant parameters. We took 500 samples from each of 
these distributions and the results are shown in Figure 4.2A alongside the relevant 
sampling distributions. 

The PoMs generated by varying the input parameters as summarized in Table 4. 1 
were solved via MATLAB stiff-solver ode15s, and allowed to equilibrate by simulat- 
ing for 500s before computing the output features. Calibrated PoMs for both the adult 
CM and hiPSC-CM were obtained by rejecting models that produced features out- 
side the experimentally observed range. We identified a collection of features from 
the literature: maximum upstroke velocity (voltage), minimum diastolic potential, 
voltage amplitude, action potential duration (APDo9), maximal departure velocity 
for Ca**, diastolic Ca?+ concentration, Ca?+ transient amplitude, Ca** transient 
time-to-peak, and time constant of Ca*+ decay. Representative scaling factors sam- 
pled from the output of the models obtained via the constraining criterion are shown 
in Figure 4.2B. In the results section, we discuss the typical features obtained and 
constraining criteria to generate calibrated PoMs for the cardiomyocytes. 
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4.2.2 Islet B-cells 


A robust secretion of insulin by the 6-cell upon glucose challenge relies heavily on 
the coupling of the metabolic oscillations in the glycolysis pathway with electrical 
oscillations (e.g. isolated action potentials and action potential bursting). Recogniz- 
ing this, the 6-cell PoM was generated via a coupling of the electrophysiological 
model developed in [9] and glycolysis model developed in [10]. A coupling of these 
models was demonstrated in [33], albeit using single mean values for the electri- 
cal and glycolytic component parameters. A simplistic illustration of this coupled 
model is shown in Figure 4.3. To delineate contributions made by the glycolytic 
and electrophysiologic systems to the overall functional heterogeneity in our 6-cell 
populations, we divide the parametric distributions into two components: ionic and 
glycolytic. 
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Fig. 4.3: Schematic representation of a 6-cell’s metabolic and electrophysiologi- 
cal circuit involved during the glucose-stimulated-insulin-secretion. Included in the 
schematic are representative formulae of the glycolytic component’s enzyme rate 
equations and a few parameters (in red) that are changed to generate PoMs. Log- 
normal distributions for two of the 14 parameters are shown in the bottom right hand 
side of the figure. 


Details of the ionic models can be found in the cited articles, but briefly electrical 
oscillations in the system were modeled as described by Equation 4.2. The maximal 
conductance (g;) of 8 ion channels were varied based on experimental data in the 
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literature (references found in [34], [35], and [36]). These variations, reported via 
coefficients of variation (op) and their sources, are summarized in Table 4.2. 

The glycolysis model comprised enzymatic reactions from glucokinase to 
glyceraldehyde-3-phosphate. The rate of ATP generation through these steps and 
via downstream oxidative phosphorylation in the mitochondria is introduced in a 
phenomenological manner via the variable, a, which represents ATP level. 


da 
VGAPDH-— kaa (4.3) 


where Vg appr is the metabolic flux due to GAPDH and k 4 is a phenomenological 
time constant.The ion channel conductance of K Arp channel is then altered via: 


g _ ÊKATP 
KATP = 
l+a 


(4.4) 


In our study, we varied the kinetics of the enzymatic reactions by varying the 
limiting rates V; max of glucokinase, phosphofructokinase, and glyceraldehyde 3-P 
dehydrogenase. In addition we varied other kinetic parameters such as the half- 
activation concentrations S wee These parameters are highlighted in red in our Figure 
4.3. Variations in these glycolytic components of our model and their sources are 
summarized in Table 4.2. This variability was again assumed to be log-normal 
distributed and implemented as such for construction of our PoMs. Sources for these 
datasets can be found in [37], [38],[39],[40]. 

Since the behaviour of individual -cells is much more diverse than those of 
cardiac myocytes, the constraining conditions are also less well-established. Thus, 
we discuss the features used for constraining the 8-cell PoMs in the results, and we 
report sub-classes of 8-cells obtained via our approach. Furthermore, we will analyze 
the effects of perturbations to the glycolysis and ionic components in maintaining 
functional heterogeneity. 


4.3 Results 
4.3.1 Cardiomyocyte POM 


As described in the methods, we allowed each POM member model to equilibrate 
for 500s before computing the output features. The last 10 action potentials of each 
simulation for these models were used to extract features that were compared against 
the constraining criteria to generate the calibrated PoMs. Constraining criteria were: 
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Table 4.2: Scaling parameters for islet 6-cell POM 
Sources to construct coefficient of variation (op) listed by the reference numbers 


Variable Mean Unit op 
Ionic component 

8&KV 1.00 nS/pF 0.37! 
&BK 0.02 nS/pA 0.80! 
2Na 0.40 nS/pF 0.49! 
Cal 0.14 nS/pF 0.34! 
&CaPQ 0.17 nS/pF 0.40! 
gcat 0.05 nS/pF 0.42! 
SKATP 0.01 nS/pF 0.897 
SHERG 0.00 nS/pF 0.303 
Glycolysis component 

VoK.max 55.60 pM/s 0.714 
Sos 8.00 mM 0.905 
VPFK,max 556.00 Mis 0.036 
SPEE 4.00 mM 0.096 
hprK 2.50 N/A 0.146 
VGAPDH,max 1.39 mM/s 0.277 


"[34], 7135], °[36], 4137], °[38], [39], 7140] 


(C1) At all times during the simulation, -100mV < Vm < 70mV, 

(C2) The change in cytosolic Ca% transient amplitude < 2% of the mean, 

(C3) The change in cytosolic Na* < 2% of the mean, 

(C4) The change in Casg transient amplitude < 2% of the mean, 

(C5) Peak Vm > —20mV, 

(C6) Action potential amplitude > 20mV, 

(C7) The standard deviation of the final 10 APDogs is less than 10% of the mean, 

(Cg) We are able to extract every feature from the model, 

(C9) The minimum diastolic potential is less than —-65mV. 
The full populations, showing both the accepted (orange) and discarded (grey) model 
configurations are shown in Figure 4.4A, together with the percent of configurations 
discarded based on each pairwise combination of constraints in Figure 4.4B. This 
figure provides a first indication of how well variability in the model behaviour is 
captured by applying a data-defined variation in the model inputs. A large number 
of configurations had to be discarded for both models, although the reasons for 
exclusion were different. Specifically, we retained only 55 of 500 models for our 
adult ventricular cardiomyocyte PoM, and 60 of 100 models for our hiPSC-CM PoM. 
For the adult ventricular PoM the vast majority of exclusions were made because the 
configuration failed either Cı (43% of configurations) or Cy (40% of configurations). 
This can also be seen in Figure 4.4A, where many configurations fail to repolarize 
below -65 mV (Cy), whereas others reach supraphysiologic membrane potentials 
(either > 70mV or < —100mV). For the hiPSC-CM model the reasons for exclusion 
were physiologically reachable but not observed in published hiPSC-CM phenotypes. 
That is, the vast majority of excluded configurations simultaneously failed criteria 
C2 - Co, each of which was chosen to address a different source of biologically 
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reachable instability. This can also be seen in Figure 4.4A, where most of the 
discarded (grey models) have attracted to a second equilibrium in membrane potential 
between —10mV and —20mV. This property of electrophysiologic bistability is a 
real biological property of cardiac myocytes [41]. It is thus a fair question to ask 
whether such configurations should be discarded for the hiPSC-CM PoM, as they 
may in fact be realistic members of most hiPSC-CM cell populations. Regardless, 
this initial characterization suggests that the experimentally-defined variation in 
hiPSC-CM input parameters results in sharp transition to gross electrophysiologic 
instabilities, most of which remain physiologically plausible. On the other hand, 
the majority of adult ventricular PoM exclusions resulted from either biologically 
unreachable phenotypes (Vm > 70mV or < —100mV), or to an unstable resting 
potential (minimum V,,, > —65mV). 


Adult PoM hiPSC-CM PoM B Adult exclusions 
50 Rais o o o oo 
Re 2200 
> 2 3200 
0 g 0 2 22030 
Wr 
-50 ka 
Aa o 2 2 Ii 
ERE 
-100 O 
0 500 1000 0 500 1000 Cı C2 C3 C4 Cs Ce C7 Ce Co 
2 hiPSC-CM exclusions 
= = 00000000 
= 1.5 CG 
= Rome © 37 37 37 36 37 37 37 37 
a rm 0 37 37 37 36 37 37 37 37 
w 1 9 
eia 0 37 37 37 36 37 37 37 37 
Y, Ries © 26 36 36 36 36 36 36 36 
0.5 Rae © 37 37 37 36 37 37 37 37 
Rites © 37 37 37 36 37 39 39 38 
0 g 37 37 37 36 37 39 39 38 
0 _ 500 1000 o 500 1000 ike © 37 37 37 36 37 38 38 38 
time (ms) time (ms) 


Cı C2 C3 C4 Cs Ce C7 Ce Co 


Fig. 4.4: A The action potential and calcium transient phenotypes from the adult CM and hiPSC- 
CM PoMs. Traces of cells accepted into the calibrated PoMs are in orange, while those discarded 
are shown in grey. B Tables summarizing the % of models excluded by the constraining criteria. 
Diagonal values show the percentage of models excluded by each criterion in isolation, whereas 
the off-diagonal entries show the % of configurations for which each combination of two excluding 
criteria were present. The meaning of each criterion, i.e. C;, is explained in the main text. 


Figure 4.5 shows the relationships between major features of the final calibrated 
PoMs for the two cardiac models, and provides a basis for comparing systematic 
differences in the phenotypic span of each PoM. The stronger relationships for the 
hiPSC-CM PoM between Ca? handling features and AP amplitude indicate the 
larger role of sarcolemmal fluxes in determining Ca? features. In contrast, the adult 
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Fig. 4.5: Scatter plots of feature pairs from both PoMs. The fitted distributions are 
shown along the diagonals using a kernel density estimator. Noteworthy differences 
between the two PoMs are boxed in the red (adult) and purple (hiPSC-CM). 


ventricular PoM Ca** dynamics are more internally determined, as demonstrated by 
stronger relationships boxed in green. These differences reflect intrinsic characteris- 
tics of the baseline hiPSC-CM and adult ventricular cardiac myocyte models. They 
are a well-known difference in the biological underpinnings of the two cell types and 
reflect the generally less “mature” Ca** handling phenotype of hiPSC-CMs relative 
to adult ventricular myocytes. 

Finally for the cardiac PoMs, Figure 4.6 shows important differences between 
the two calibrated PoMs in terms of degree to which the variability in their AP 
and Ca? handling features reflect those reported in the literature. Specifically, it 
is clear that the adult ventricular PoM features (top row, orange histograms) are 
generally more variable than is reported experimentally (dashed lines). Only the 
APDop feature for this POM appears to approximate reported experimental variation. 
In contrast, the hiPSC-CM feature distributions (bottom row, orange histograms) 
are less variable than published reports of these phenotypes (dashed lines). This 
presents a fundamental distinction, and suggests that literature reports of variation 
in the underlying currents and fluxes overestimate the true physiologic variability in 
human adult cardiac ventricular myocytes. In contrast, the same techniques applied 
to hiPSC-CMs underestimate observable phenotypic variability. We explore this 
important observation in the discussion. 
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Fig. 4.6: Histograms showing the variability in the features of the PoM. The first 
row shows the adult ventricular cardiomyocyte PoM, and the second row shows the 
hiPSC-CM PoM. In each plot, the probability density function of a representative 
Gaussian distribution shows the variation typically found in measurements of these 
features in the literature. These distributions have been normalised such that they are 
centered on the means of each histogram. 


4.3.2 Islet B-cell POM 


The individual models for the 6-cell PoM were simulated to capture cell behaviour 
over 10 minutes - the acute phase of glucose-stimulated insulin release. The param- 
eter variations provided in Table 4.2 were used to generate 1000 cells to assess the 
importance of 3 sources of heterogeneity, by varying model parameters associated 
with: (a) only glycolytic (enzyme) model components, (b) only ionic (ion channel 
and transporter) components, (c) both ionic and glycolytic components. In all three 
cases, each model configuration was simulated under high glucose (10 mM) and low 
glucose (2 mM) conditions. Thus, the total number of simulations across all 3 cases 
(3000 cell model configurations), at both glucose concentrations, was 6000. 

Figure 4.7A shows V,, and cytosolic calcium for representative simulations at 
high glucose. We see cells 1 through 3 demonstrate similar features as experimentally 
observable sub-classes of 6-cells, i.e. spiking, bursting, and plateau phenotypes. Non- 
firing cells such as those shown in cells 4 and 5 were also observed. As for the cardiac 
PoMs, in order to generate a calibrated 6-cell PoM, we have to invoke and apply 
a set of constraining criteria to reject configurations that produce features outside 
the ranges of experimental observations. However, because the 6-cell phenotype 
is the most heterogeneous of all cell types assessed here, there exists far fewer 
such objective constraints in the 6-cell literature. However, some benchmarks are 
available. For example, [42] suggested two metrics of activity based on B-cell V, 
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Fig. 4.7: Representative simulation outputs and categorizations. A. Representative 
simulation outputs of membrane potential and cytosolic calcium representing dif- 
ferent sub-classes of B-cells. 1. Spiking cells 2. Bursting cells 3. Plateau cells 4 & 
5. Non-firing cells. Scales on the plots represents time of 240s. B. Plot of activity 
fraction vs. Apeak for the 5 representative cells. C. Number of peaks in the Fast 
Fourier Transform (FFT) of the membrane potential traces for the 5 representative 
cells. 


(the activity fraction and Apeak), as means of classifying phenotypic sub-classes. 
Activity fraction was defined as the fraction of time V,m was above an arbitrary 
threshold level. The definition of Apeak was set as the difference of the two major 
components of V,,,-time probability distribution (i.e. the difference between “resting” 
Vm and “active” Vm). As for [42], we did not observe clear delineation of the three 
sub-classes, nor could we readily differentiate firing vs. non-firing -cells via these 
metrics (Figure 4.7B). However, we surmised that Fast Fourier Transform (FFT) of 
the Vm signal may allow more clear distinction of these classes. Figure 4.7C shows 
the number of FFT peaks corresponded with the number of dominant harmonics in 
the V,,, signal for each representative recording in Figure 4.7A. 

Based on this analysis, we set the constraining criteria as follows. Any cell with 
less than 100 peaks in the Vm peaks, was considered to be quiescent and thus 
discarded. For cells that had more than 100 peaks in their voltage signal, FFT was 
performed on the voltage signal to classify them into spiking, bursting, or plateau 
cells. Based on our initial analyses, cells with less than 20,000 peaks in the voltage 
FFT were classified as plateau cells. Bursting cells were classified as those having 
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anywhere above 20,000 FFT peaks but less than 40,000 peaks. Finally, cells with 
FFT peaks greater than 40,000 were classified as spiking cells. 


om 
iv] 


viable cells) 


40.5 


% Viable cells 


B- subtype (% 


28.7 0.0 
| 19 15.5 a 200 
a E Li 0 i G LG HG LG 


2 Component lonic Glycolysis 


Cc D 
0.4 0.4 <0.0001 
<0.0001 


0.1569 0.0012 200000: 
0.0142 <0.0001 rr 
ü < EE 150000: 
= 
a 
kal 
3 © 100000. 
0.2 i : 
Fi 
8 
50000: 


n=322 n=165 n=179 


HG LG 
2 Component lonic Glycolysis 


<0.0001 


<0,0001 0.0001 


o 


Mean cytosolic Ca?* (uM) 
Peak Cytosolic Ca?*(uM} 


° 
Peaks (Cyt 


B 


n=322 n=165 n=179 


n=322 n=165 n=179 


0.0 
Spikers Bursters Plateau Spikers Bursters Plateau Spikers Bursters Plateau 


Fig. 4.8: A. Percentage of viable cells based on the chosen constraining conditions 
for six different cases (1000 runs each). B. Distribution of the sub-classes of B-cell 
from the viable candidates for each of the six cases. HG: High glucose (10 mM). 
LG: Low glucose (2 mM). Analyses of the cytosolic calcium traces of each of the 
sub-classes of 6-cells in the viable population obtained with high glucose: C. mean 
values, D. peak values, E. Number of FFT peaks. Statistical differences between the 
groups were determined by one-way ANOVA followed by Tukey’s HSD post-test. 


Based on these constraining criteria, we quantified the viable cells that were not 
rejected for both glucose concentrations and all 3 cases of heterogeneity (glycolytic, 
ionic, both), as shown in Figure 4.8A. The greatest reduction in viable cells was 
observed for glycolytic heterogeneity only, thus highlighting the narrow metabolic 
parameter range functional 6-cells operate under and its importance for achieving 
robust insulin secretion. Another interesting aspect of these data was the number 
of active cells at low (“resting”) glucose (2 mM). This feature of basal activity 
is unique to human f-cells, and underlies basal insulin release at the lower blood 
glucose concentrations present in humans. Thus, our PoM approach was able to 
capture this aspect of human-specific B-cell physiology. 

Further analysis of the viable candidates allowed us to assess the fraction of each 
of the 3 PoMs belonging to the sub-classes of 6-cell phenotypes at each glucose 
concentration (Figure 4.8B). The lowest fraction of bursting cells was observed 
in the low-glucose (LG) cases when both components or only ionic components 
were varied. Upon altering glycolytic components alone, plateau cells comprised 
the only viable sub-class, thus suggesting the importance of coupling the ionic and 


4 Functional Heterogeneity 59 


glycolytic elements to generate the full range of functional phenotypes observable 
in human £-cell populations. Mean cytosolic calcium of the viable cells within 
our PoM were not statistically different between spiking and bursting cells, sugges- 
tive of less pronounced basal shift of calcium signalling (Figure 4.8C). However, 
plateau cells exhibited a mean cytosolic calcium signature distinct from the other 
two classes, while peak cytosolic calcium grouped plateau and bursting phenotypes. 
(Figure 4.8D). Applying FFT analysis to the calcium signals also showed that the 
number of harmonics in the signal could distinguish the sub-classes of B-cells, thus 
suggesting that this may be a viable means of distinguishing functional heterogeneity 
in experiments, for which calcium recordings are considerably easier. 


4.4 Discussion 


In this manuscript, we report the construction of PoMs for three electrophysiologi- 
cally active tissues that were calibrated against experimental data from the literature. 
PoMs have proven useful in investigation of cardiac electrophysiological variability 
and recent studies have furthered the methodology by explicitly incorporating exper- 
imental data into the construction of populations of models [43, 6]. In our CM PoMs, 
we observed that the fluctuations in the input parameters for the models allowed us 
to generate cells with features comparable to those observed in both primary and 
hiPSC-derived CMs, albeit with differing coverage of the range of function observ- 
able in experiments. In both cardiac PoMs we assumed the main source of variability 
is differing expression of ion transporters, and that this variation is responsible for the 
observable variation in their Vm and calcium cycling features. In hiPSC-CMs par- 
ticularly, this variation is likely to be pronounced given their immature (“fetal-like”) 
phenotype and the degree to which this can be matured. Importantly, we observed 
that experimentally-defined variation in the input parameters of the cardiac models 
resulted in functional output that differed in its range for the two models. That is, 
for the adult CM PoM, the range of AP and calcium handling features was broader 
than observed experimentally, whereas it was narrower than experimental for the 
hiPSC-CM PoM. This is a fundamental observation about the agreement between 
two types of data recorded in these cell types. For the adult CM it could be that bio- 
logical covariance (which we have not applied in PoM construction) among the ionic 
currents reduces the range of reachable AP and calcium features in real cells, or that 
the relative scarcity of adult human CM data means that these two classes of data are 
not yet fully internally consistent. In contrast, the greater experimental variability 
for hiPSC-CM phenotype features suggests either that experimental ionic current 
measurements are subjected to overly constraining inclusion criteria. Alternatively, 
our hiPSC-CM model formalism may somehow be constrained to be overly stable 
with respect to the phenotypic variation resulting from variation in the input currents. 
Regardless, these questions are fundamental to understanding how perturbations to 
either adult or hiPSC CMs should be expected to impact those tissues in vivo. 
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The islet -cell PoMs were generated via alterations in both the parameters in- 
volved in glycolysis model components as well as the ionic conductance components. 
We observed that the number of harmonics in V,, of the simulated computational 
cells was enough to distinguish different sub-classes of B-cells observed experi- 
mentally, although we believe the proposed constraining criteria requires further 
validation via prospective screening of experimentally obtained data sets. Further- 
more, we observed that the functional heterogeneity of the 6-cells is dependent 
on the coupling of the metabolic and ionic components, as noted by the number 
of viable cells and their sub-types when only glycolytic components were varied. 
This, in part, explains why glycolytic bottlenecks in stem cell-derived -cells show 
non-robust glucose-stimulated-insulin-secretion behaviour [44, 45]. The calibrated 
PoMs generated in this study were able to capture electrical oscillations at low glu- 
cose, which is a unique feature of human £-cells (compared to rodents). Thus, the 
methodology developed should be applicable for creating human-specific 3D models 
of coupled f-cell clusters with and without partner cells such as a-, 6-, €- and pan- 
creatic polypeptide cells. When coupled with experimental studies, these 3D models 
could be invaluable to understand the effect of heterogeneity on intact islet function, 
and to develop sophisticated, robust tissue-engineered platforms for therapeutics. 

Here we have demonstrated the use of experimental constraints to first construct 
and then calibrate PoMs of three human cell types. The cardiomyocyte PoMs partially 
reproduced observable variation in the functional features for their respective cell 
types, albeit with clear systematic inaccuracies. The ability of the human $-cell POM 
to recapitulate the range of function observed in experimental 6-cell recordings 
is encouraging for applying this strategy to integrated function of B-cells in human 
islets. Constructing new methods for constraining these PoM approaches will require 
continued advancement of single cell assays. 
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Abstract The tripartite synapse or “neural threesome” refers to the interplay in 
the synapse between neighbouring neurons, the synaptic cleft, and the surrounding 
glial cells. Despite extensive research, the effects of glial cells, such as astrocytes, on 
signal transduction between neurons are not fully understood. The Kirchhoff-Nernst- 
Planck (KNP) and Extracellular-Membrane-Intracellular (EMI) models constitute a 
promising framework for modeling these kinds of systems. However, they lack the 
neurotransmitter-related mechanisms that are necessary to bridge signal transduction 
across the synaptic cleft. Here, we propose an extension to the KNP-EMI model by a 
spatio-temporal diffusion-based description of the most prominent neurotransmitter, 
glutamate, that allows for investigation of the contribution of astrocytes to the func- 
tionality of the synapse. We validate our model by showing that the presence of an 
astrocyte in the domain affects the glutamate flux across the postsynaptic terminal, as 
observed physiologically. The proposed extension offers a sufficiently simple way of 
integrating synaptic glutamate dynamics into the KNP-EMI framework. It introduces 
the relevant interactions between electrical activity and diffusion processes at the 
tripartite synapse that are necessary to assess how astrocytes might contribute to the 
functionality of the synapse. This work has implications for future studies involving 
glial mechanisms and other charged species within the KNP-EMI framework. 
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5.1 Introduction 


Circuits in the brain are composed of complex connections of both neurons and 
glia. Neurons communicate via patterns of electrical signals that are propagated by 
molecules and ions. These signals, called action potentials, communicate information 
that gives rise to cognitive function, behaviors, and movements. Action potentials 
travel from neuron to neuron across small gaps called synapses. The neuron sending 
the signal, the presynaptic neuron, releases specific molecules called neurotransmit- 
ters to the receiving neuron, or the postsynaptic neuron. However, neurons are not 
the only cells that are found in the synapse. A type of glial cell called an astrocyte 
envelopes the synapse and alters information transmission by modulating the neuro- 
transmitters as they diffuse across the synapse [1]. Thus, together, the neurons, the 
synapse, and the astrocyte form a “neural threesome” [2], or tripartite synapse. The 
role of this astrocyte mediated modulation in electrical and chemical interplay is not 
well understood. 

One of the ways astrocytes are thought to affect synaptic transmission is through 
reuptake of the neurotransmitter glutamate. Once an electrical signal reaches the end 
of the presynaptic neuron, prepackaged vesicles release glutamate into the synaptic 
cleft. After release, glutamate diffuses across the synapse where it binds AMPA 
receptors on the postsynaptic neuron, resulting in a net inward sodium (Na*) current 
through the ligand-gated AMPA receptor ion channels. As a result, the intracellular 
potential starts to increase which triggers the opening of nearby voltage-gated Na* 
channels. This channel opening allows for a larger Na* influx, bringing the postsy- 
naptic neuron closer to its threshold for depolarization. If the signal is strong enough, 
subsequent ion dynamics across the postsynaptic membranes result in the onset of a 
new action potential and signal propagation continues. 

As glutamate diffuses across the synapse, astrocytes take up a percentage of the 
glutamate ions through high-affinity glutamate transporters [3, 4]. This mechanism 
is suggested to be a way of reducing cellular toxicity, as too much glutamate in the 
synapse can have excitotoxic effects [5]. Recently, it has been found that astrocytes 
also release glutamate to neurons in response to extracellular glutamate cues [6]. 
Because of this dual mechanism of glutamate uptake and release, astrocytes act as a 
buffer to balance excitatory and inhibitory neuronal activity. 

Here, we propose a simple mechanistic model of glutamate in the tripartite 
synapse that allows for the assessment of these questions within the Kirchhoff-Nernst- 
Planck (KNP) - Extracellular-Membrane-Intracellular (EMI) modeling framework. 
The KNP-EMI model is a mechanistic mathematical model for excitable tissue in 
spatiotemporal resolution. It is stated in the form of a set of coupled partial differen- 
tial equations, and includes interactions between ionic transmembrane currents and 
electric potentials [7, 8, 9]. In [7], Ellingsrud et al. show that their KNP-EMI model 
accurately describes depolarization along excitable membranes. Thus, it constitutes 
a truly spatio-temporal model of signal transduction along a cell, for example, along 
an axon or a dendrite of a neuron. However, the model in its current form does not 
account for neurotransmitters in the extracellular space, nor for neurotransmitter- 
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related mechanisms and interactions. For this reason, the KNP-EMI model is not 
able to capture the regulation of signal transduction at the synapse. 

We propose an extension by introducing glutamate to the model. In order to bridge 
the synaptic gap, we identify three relevant biophysical mechanisms: glutamate 
release into the cleft by exocytosis from the presynaptic neuron, diffusion in the 
extracellular space, and binding to postsynaptic AMPA receptors. All three of these 
processes are considered in our implementation. In addition to this, we assume 
that glutamate uptake is the fundamental interaction of the astrocyte at the synapse 
whereby astrocytes regulate signal transmission. Therefore, our proposed glutamate 
diffusion model describes the properties of signal transduction unique to the synapse; 
in particular, how astrocytes affect the glutamate gradient and how a changing 
glutamate gradient affects signal propagation. Our extension to the KNP-EMI model 
features a simplified yet qualitatively complete description of signal transmission 
across neurons. The extended model will allow for a comprehensive investigation of 
the astrocyte’s role in signal propagation across neurons in a way that takes spatial 
effects, such as explicit geometries and local concentrations, into account. 

The EMI model constitutes a general framework for modeling excitable tissue 
ranging from cardiac cells to neural systems [10]. Besides its broad range of appli- 
cation, it allows for the incorporation of explicit morphologies within the system of 
interest. In combination with the KNP framework, the KNP-EMI model allows for 
the representation of a complicated geometry and detailed electrochemical profiles 
[8]. In order to leverage its power, we design our glutamate model under consid- 
eration of its integrability into the KNP-EMI model. Although its formulation is 
independent of the respective morphology, our implementation assumes a simplistic 
geometry in order to show its functionality as a proof of concept. 

Our paper is structured as follows. First, we familiarize the reader with the ge- 
ometry of the problem conceptually and derive the glutamate model. Subsequently, 
we demonstrate a way to estimate the model parameters based on literature val- 
ues, restricting our model to hippocampal synapses. After providing implementation 
details, we validate the glutamate model by showing that it behaves as intended. Fi- 
nally, we discuss the implications of our model within the context of the KNP-EMI 
framework. 


5.2 Methods 


The glutamate model consists of two parts. The first part describes the transport of 
glutamate across the synaptic cleft by diffusion, whereas the second part describes 
how AMPA receptors react in response to glutamate to trigger a Na* current across 
the postsynaptic terminal. Here, we describe the computational domain and the 
governing equations for the glutamate model. We then describe the parameters that 
are used within the model, and finally, the numerical methods used in our simulations. 
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5.2.1 Representation of the Computational Domain 


The tripartite synapse features a variety of membranes that behave differently. The 
characteristic behavior of each membrane is determined by the density and type of 
functional components it exhibits, e.g. channels, pumps, and receptors. Figure 5.1 
shows a schematic of our envisioned computational domain for the full KNP-EMI 
model. It resembles a simplified and abstracted two-dimensional realization of a 
tripartite synapse setup. The domain is divided into two distinct subdomains: the 


Pastrocyte 

Thresynaptic neuron 
Thostsynaptic neuron 

I, (not glutamate permeable) 


Open boundary 
Glutamate transporter 


Hodgkin Huxley ion 
channel 


AMPA receptor 
Intracellular space 


Extracellular space 


Fig. 5.1: Representation of the computational domain in our glutamate KNP-EMI 
model. Synthesis of computational domains, channel dynamics, and boundary con- 
ditions. 


extracellular space, Qe, and the intracellular space, Qi. Qij represents the interior 
of the postsynaptic neuron’s dendrite. The interface boundary separating Qe from 
Q; consists of the postsynaptic terminal, Ipost; and extrasynaptic membrane of the 
dendrite, I4. We assume that I4 is impermeable to glutamate. In contrast, Ipost is a 
region with an affinity for glutamate due to the presence of AMPA receptors, which 
are only located on the part of the postsynaptic membrane that faces the synapse [11]. 
The exterior boundary consists of three biophysically distinct types. pr. denotes 
the membrane between the presynaptic neuron and the synapse. I, represents the 
membrane on the astrocyte, which also has an affinity for glutamate due to the 
presence of glutamate transporters that clear it from the synapse [4]. The remaining 
part of the exterior boundary is open for glutamate to pass, since the neurotransmitter 
freely diffuses through the extracellular space out of the synaptic region. We note 
that because this report only covers the glutamate model, our computational domain 
is for now restricted to the extracellular space, Qe. 
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5.2.2 Mathematical Modeling of Glutamate Dynamics 


The glutamate transport is modelled using a pure diffusion equation: find the gluta- 
mate concentration g = g(x,t) such that 


g=DeAg on Qx[0,T] . (5.1) 


Here, De is the diffusion coefficient for glutamate within the extracellular space. 
The glutamate release from the presynaptic terminal and absorption of glutamate on 
the various cell membranes are modelled using the following combination of von 
Neumann and Robin boundary conditions: 


-D.Vg-nr, = 0 on Tax [0,7] (5.2) 
-D.Vg-nr, =  Kuptake glr, on I,x [0,7] (5.3) 
-De V8 Mtg = Kvind Slr sox on [post X [0, T] (5.4) 
—De Vg nr =— Diy 6(t-t)) on Tprex[0, 7]. (5.5) 


Note that here, ng denotes the outward pointing normal vectors of the respective 
boundary surfaces. As a consequence, the left hand sides denote the outward fluxes of 
glutamate through them. We model glutamate binding to the AMPA receptors as well 
as glutamate uptake at the astrocyte’s membrane by first order reaction kinetics. Thus, 
Kvina Specifies the affinity of the AMPA receptor towards glutamate and Kuptake the 
absorption rate. The reflecting boundary condition at I'y in Equation (5.2) realizes the 
impermeability of glutamate through the extrasynaptic membranes of the dendrite. 
Equation (5.5) describes the source of glutamate in our model, which is vesicular 
release due to exocytosis. In accordance with Clements et al. [12] and Jonas et al. 
[13], we model an instantaneous release of vesicles at release times f;. vi, denotes the 
vesicular spatial profile defined in Equation (5.9). On the remaining open boundaries, 
we assume homogeneous Dirichlet boundary conditions 


aire 50 ¢ (5.6) 


This choice is made for simplicity, and future work will benefit from more realistic 
modelling of free diffusion across the open boundaries. Since glutamate is supposed 
to be cleared from the synaptic cleft after exocytosis, we assume an initial equilibrium 
of g(x,t =0) =0. 


5.2.3 Modeling AMPA Gating Dynamics 


We model the dynamics of the AMPA receptors on the postsynaptic terminal using 
a model proposed by Tewari and Majurmdar [14]. The ODE model describing the 
AMPA gating dynamics reads 


70 EMI for Astrocyte-Neuron Interactions 


d 


J; AMPA = = Campa 8|r (1 -mampPa) -PAMPA MAMPA - (5.7) 


Here, the glutamate concentration, g, couples linearly to the opening rate of the 
AMPA receptor, «ampa. The ODE model governs the dynamics of the gating vari- 
able, mampa, point-wise along the boundary Tpost, and depends on the local glutamate 
concentration, Bit . We intend to initialize the system in its equilibrium state. The 
initial state g(x, t= 0) = 0 Cortasponds to the equilibrium state of the glutamate 
model, For g =0, we get m,;p, =0 as the ae state of the gating variable 
MAMP a ât t = 0. Therefore, we choose m AMPA(f = 0) = MANIP a =0. 

The gating variable, mampa, controls an AMPA iS Na* current, TAMPA, 

through the postsynaptic terminal. This current is given by 


TAMPA = ZAMPA MAMPA (Vpost— Vampa) on Tpost (5.8) 


[15, 14], with Vpost being the local membrane potential at the postsynaptic neuron 
and Vampa referring to the reversal potential. In the KNP-EMI model, Vampa is 
associated with the Na*-specific Nernst potential. 


5.2.4 Estimation of Model Parameters 


All parameters for the AMPA receptors were chosen in the same way as the original 
model by Tewari et al. [14]. Following Clements et al. and Tewari et al., we assume 
Las with a diameter dyes = 40nm, containing glutamate at concentrations of 
cS! = 60mM. In addition, we assume a homogeneous glutamate surface density 


A across the vesicle’s cross section on the presynaptic terminal. We may now 


define a vesicular profile at Ipre, vi, as 


(S. V [y — Yves| < dves 


5.9 
0 else, On 


vi (y) = 


where Yvyes denotes the center position of the vesicle on the membrane and y denotes 
the coordinate parameterizing the one-dimensional boundary surface. We assume 
that the arrival of the action potential at the presynaptic bouton at any time t = trel 
results in the release of one vesicle. On release, the total mass of glutamate contained 
within one vesicle exocytosed across the membrane [ye within the period [że] — 
€, trey + €], € > 0 can be written as 


felt € 
M= a f De Vg- nr, dS 


trel— € Thre 


Yvestdves /2 
= J ves | dy = one dyes 
Yves — Aves /2 
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. ; 2 GI _4,p3 „Gl . . 
On average, this total mass is Myes = Vves Cyes = 3 T Rig CVey assuming spherical 


vesicles with radius Ryes = dyes/2 . This yields a mean surface density of 
m 
Ves = g Cves Aes - (5.10) 
Released glutamate diffuses throughout the synaptic cleft. In accordance with 
[13], we choose a diffusion coefficient of De = 3-10°nm?ms~! for glutamate in 
extracellular space. In common reaction kinetic models of synaptic glutamate, the 
total clearance rate amounts to ke = 10ms™! [15, p.4]. 4/5 of synaptic glutamate is 
cleared by glial uptake, which leaves only 1/5 being removed by processes on the 
postsynaptic terminal including uptake and binding to receptors [16]. Furthermore, 
glutamate gets cleared passively by diffusing out of the cleft. We use these propor- 
tions to estimate the surface uptake and binding rate densities Kuptake, Kvina. Since 
De >> ke, diffusion is the dominant process governing the synaptic glutamate dis- 
tributions. Thus, we assume De — œ locally when focusing on the reaction kinetics 
on the membranes. In this limit, the system is well-mixed and glutamate is homo- 
geneously distributed such that g = go € R. It behaves like a single-compartment 
model, which allows us to directly compare the respective flux terms. In particular, 
we require that all of the clearing fluxes add up to kc, by 


1 
5 ke 8c = ef Kvina dS (5.11) 
post 
4 
5 ke 8c = 8c Kuptake dS . (5.12) 
Ta 


We assume that the membrane properties do not change locally, therefore, we assume 
homogeneous surface rate densities. Thus, we have simple expressions for both 
surface rate functions, 


(5.13) 
(5.14) 


Here, ws denotes the width of the synaptic cleft and h; describes the full height of the 
synaptic terminals. Table 5.1 summarizes all parameters and their values associated 
with our glutamate model. 


5.2.5 Numerical Methods 


In our validation studies, we solve the glutamate diffusion model numerically using 
the finite element method with piecewise continuous linear elements for spatial 
discretization and an implicit Euler scheme for temporal integration. 
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Table 5.1: Parameters and their values relevant for the glutamate model 


Symbol Value Unit Description 

De 3-105 nm? ms™! Diffusion coefficient in the extracellular space 

Kvind 5.107? nm`!ms™! AMPAR specific glutamate binding surface rate density 
Kuptake 4-107! nm! ms™! Glutamate uptake surface rate density of the astrocyte 
i. 5-105 mM nm? Glutamate release surface density 

dy 40 nm Mean diameter of a vesicle 

Ws 20% nm Width of the synaptic cleft 

hy 400 “ nm Full height of the synaptic terminals 


a Values according to Clements et al. [12]. 


Let Vp denote the set of continuous piecewise linear functions restricted to a 
discrete triangulation of Qe, such that v € Vp vanish at the open boundaries Vopen. 
Then, the discrete variational form of the glutamate diffusion model at time t = n At, 
where n = 0,1,..., and Ar is the time step, reads: find g € Vp such that 


1 
/ E8 +e Vvde | Kkuptake glr, vass | Kpina 8lriou V dS 
Qe Ta post 


_ 1 0 i ok 
-f ne vaD fobs (t-t;)vds Wve Vp, 


Here, g? refers to the solution at the previous timestep tọ = (n — 1)At, and 6* refers 
to a distribution that, similar to the dirac delta in Equation (5.5), has the property 
that the glutamate released from each vesicle coincides with its total mass when 
integrated in time. This distribution may, for example, be a box function which 
evaluates to 1/6, over some interval [t,t+6,], or a Gaussian function parameterized 
by the length scale 6;. 

On the other hand, we implement the AMPA model (5.7) using an explicit Euler 
integration scheme: 


(5.15) 


0 0 0 0 
MAMPA = M AMpa + At (campag lia (1 = mima) — PAMPA Mara) . (5.16) 


The model is implemented and solved using FEniCS [17]. 


5.3 Results 


In the following section, we present the results of our validation study regarding 
our glutamate model. Figure 5.2 illustrates our experimental setup. It shows the 
glutamate concentrations at three different time points of our numerical solution 
in the synapse. We start the simulation with the release of a vesicle from Tyre at 
t = 0.01 us. In our study, we compare two cases: glutamate release located close to 
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the site of the astrocyte at the bottom of the presynaptic terminal, and release at the 
opposite site of the terminal, far from the astrocyte. Figure 5.2 displays both cases. 
At t = 0.01 us, upon release of the glutamate into the synaptic cleft, Figure 5.2 
shows highly localized concentrations of glutamate near the regions of release at 
the presynaptic terminal, indicated by the red regions of high contrast. After 0.1 us, 
the pronounced glutamate concentration profiles have become smoother and the 
distribution has significantly broadened. At t = 0.5ys, the glutamate distribution has 
become diffusive, and it is clearly visible that glutamate is cleared from the synaptic 
cleft. It should be noted that the depiction in Figure 5.2 is compressed in height. In 
fact, the synapse is higher than it is wide by a factor of 20. Note that the apparent 
increased transport along the horizontal axis is an artifact of this compressed domain. 
In contrast to our expectations, there seems to be no apparent difference between 
0.01 0.1 
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21 
18 
14 
-10 
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Fig. 5.2: Glutamate dynamics in the synaptic cleft. Glutamate distribution in the 
synaptic cleft is shown at three different time points. Vesicles are released at t = 
0.01 us on the presynaptic terminal (left boundary) either near the astrocyte’s location 
(bottom) or far away from it near the open boundary (top). Glutamate bridges the 
postsynaptic terminal (to the right) by diffusion. For the sake of visualization, the 
synapse is compressed in height. 


Time t [ys] 


Glutamate concentration [mM] 


the two distributions suggesting that the astrocyte’s influence on the clearance of 
glutamate might be negligible compared to diffusion. 

Figure 5.3 summarizes the glutamate fluxes through the membranes of interest 
and the resulting AMPA-related Na* current in a quasi-static approximation of the 
postsynaptic membrane potential. Figure 5.3A shows the glutamate flux J Astro across 
the astrocyte’s membrane T; as a function of time. After the release near the astrocyte 
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(golden curve), the net flux of glutamate across the membrane increases rapidly and 
starts to decrease again after 0.25 us as the glutamate gets cleared from the synaptic 
cleft. In contrast to that, we do not measure a significant glutamate uptake through 
the astrocyte at all if the vesicle is released at the opposite end (teal curve). We note 
that shortly after the release, Jastro becomes negative. According to the boundary 
condition in Equation (5.3), this corresponds to non-physical negative glutamate 
concentrations at the boundary. Since g > 0, this observation indicates numerical 
instabilities that may be attributed to the Gibbs phenomenon. Despite that, the 
concentration and Jastro quickly stabilize and exhibit the expected behaviour. 
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Fig. 5.3: Flux dynamics vs. time in the glutamate model. a): Total glutamate uptake 
by the astrocyte. b): Total glutamate flux through the postsynaptic terminal, Jpost, due 
to binding to the AMPA receptor. c): Total Na* current Zappa as a result of AMPA 
activation caused by the arriving glutamate in case of a quasi-static membrane 
potential at Vpost = —65 mV. 


Figure 5.3B depicts the glutamate flux through the postsynaptic membrane, JPost» 
in time. JPost has important implications for the excitation of the postsynaptic neuron 
because it indirectly reflects the activation of the AMPA receptors that control Na* 
channels on the membrane. The profile of glutamate uptake through the postsynaptic 
terminal is qualitatively similar to that of the astrocyte. On the falling branch of 
the curve however, we see that if the astrocyte is locally present, less glutamate 
reaches the postsynaptic terminal as a result of active clearance by the astrocyte. 
This validates our model with respect to the literature [3, 5, 4]. 

Figure 5.3C shows the magnitude of the AMPA-regulated Na* current across 
the postsynaptic terminal in time. Thus, it reflects the total activity of the AMPA 
receptors on [post as glutamate approaches that membrane and binds to the receptors. 
Tampa is calculated from the AMPA gating variable, mampa, using Equation (5.8). 
We assume that in the simulated period, the postsynaptic membrane potential has 
not deviated too far from its resting potential despite the finite AMPA current. With 
this quasi-stationary approximation, we justify Vpost ~ —65mV, the neuron resting 
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potential, and are able to give a coarse estimate of Zampa Without having to couple 
our model. Similar to Tewari et al. [14], we choose gampa = 0.35nS and Vampa = 
OmV. Negative signs indicate a net Na* current out of the dendrite instantiating 
depolarization across the membrane. As glutamate reaches the postsynaptic neuron, 
we see a growing increase in the AMPA-specific Na* current, as expected. In contrast 
to JPost Zampa does not saturate within the simulated time window. We note that 
while less glutamate reaches the postsynaptic neuron in the region of influence of 
the astrocyte (Figure 5.3B), unexpectedly, the resulting AMPA current Zampa of this 
region appears to be more pronounced. 


5.4 Discussion 


We developed a diffusion-based, spatio-temporal model of the neurotransmitter glu- 
tamate in glutamatergic hippocampal synapses. We carried out validation studies 
demonstrating that our model reproduces the astrocyte’s role in glutamate uptake 
in the tripartite synapse. In addition, we observed that astrocytic glutamate uptake 
translates into a decrease in glutamate flux across the postsynaptic terminal, show- 
ing that astrocytes temper the excitatory signal that is transmitted to postsynaptic 
neurons. This is a widely recognized key role of astrocytes [3, 5, 4]. We also saw 
in our model that the binding of glutamate to AMPA receptors in the postsynaptic 
neuron results in a net AMPA-specific Na* current Zampa, initiating depolarization in 
the postsynaptic neuron. Thus, our model successfully translates synaptic chemical 
cues into an electrical signal in the postsynaptic neuron. Unexpectedly, we observed 
that the total postsynaptic glutamate flux, Jpost, inversely affects the strength of the 
AMPA-specific Na* current. This effect could be explained by the astrocyte causing 
a narrower glutamate distribution, resulting in a locally restricted but higher AMPA 
activation. As a consequence, the total Na* current might be stronger. The glutamate 
profiles in Figure 5.2, however, indicate otherwise. Moreover, our quasi-static ap- 
proximation of Zampa does not resemble the actual Na* current dynamics accurately. 
Using the correct dynamic potentials, e.g. by coupling our glutamate model to the 
KNP-EMI model, might yield different results. Finally, we observed characteristics 
in our simulation results that indicate issues of numerical stability near the astrocyte 
during release. Similar instabilities might be causing the discrepancy in Zampa that 
we can not explain in the scope of our model. Thus, the next logical step would be 
to study the stability and convergence of our methods comprehensively. 
Furthermore, we observed a significant difference in postsynaptic activity with 
and without the astrocyte. However, the effect is not as pronounced as we would 
have expected. Conversely, it appears rather minute. The reason for that might lie 
in our choice of the astrocyte’s location in terms of distance from the neurons. 
In our model, we have chosen a distance of 50 nm, which is 2.5-fold larger than 
the width of the synaptic cleft. As a next step, we suggest situating the astrocyte 
closer to the neurons. Lastly, the astrocyte’s role of clearing the synaptic cleft is 
expected to be most significant when it is saturated with glutamate, for example, 


76 EMI for Astrocyte-Neuron Interactions 


in bursting neurons. In this case, excessive concentrations of neurotransmitters can 
have a neurotoxic effect. In our study, we only considered single-vesicle release, but 
our model accounts for an arbitrary amount. Consequently, with longer simulated 
time windows, multi-vesicular release studies resembling spike trains of arbitrary 
frequency are an obvious and promising connecting factor for future work. 

Other computational studies measuring glutamate flux in a synapse could not be 
found in the literature. However, some physiological studies may be used to con- 
textualize our results. Although glutamate flux across an astrocyte is not directly 
measured in vitro, it can be carried out indirectly by measuring the changing glu- 
tamate concentration of a medium containing astrocyte cultures. In this regard, the 
studies by Farinelli and Nicklas [18] and Fonseca et al. [19] show similar glutamate 
clearance dynamics as our results (e.g., Figures 5.2 and 5.3, respectively). These 
astrocytic clearance dynamics are measured over a much longer time period. How- 
ever, the observed pattern may be extended to our study in that there is an initial 
fast rate of glutamate clearance, followed by a plateau stage during which the rate of 
clearance slows as the relative glutamate concentration approaches 0 mM. 

A clear extension to our study would be utilizing the model in tandem with 
the KNP-EMI model. Our model design integrates well into EMI or KNP-EMI 
frameworks by adding the AMPA-specific contribution to the Na* current, ZAMPA. 
In the EMI and KNP-EMI frameworks (e.g. [8]), the choice of the functional form 
of Jion is subject to modelling, usually in form of a solution to a Hodgkin-Huxley 
type ODE system [20]. When integrating our model into the KNP-EMI framework, 
we suggest setting Jion = Zampa On the postsynaptic terminal (Tpost in Figure 5.1) and 
Hodgkin-Huxley based models on the remaining postsynaptic membranes (Ta). 

Considering further astrocyte-related mechanisms, such as the dual function of 
glutamate uptake and release [6], the KNP-EMI model is expected to be especially 
well-suited for modelling glutamate uptake. As a first step into modeling neuro- 
transmitters within the KNP-EMI framework, we designed our model with simplicity 
in mind. For this reason, our model only features what we identified as the most 
relevant processes for signal transmission, excluding presynaptic glutamate re-uptake 
and astrocytic glutamate release. These mechanisms should, however, be considered 
in a refined model, because they are supposedly relevant for realizing glutamate 
buffers in the synapse. Additionally, the dual function mechanism is controlled 
in part by intracellular pH levels [21], which may be well handled by the ionic 
concentration and diffusion components of the KNP-EMI model. 

In summary, we successfully developed a glutamate extension to the KNP-EMI 
model in the tripartite synapse. Our model exhibits all the relevant mechanisms 
necessary for signal transmission between neurons across the synaptic cleft. In our 
model, the astrocyte modulates the signaling molecule concentrations as expected. 
Their effect could be more pronounced, but we expect that our choice regarding the 
astrocyte’s distance from the synapse had a significant influence on this outcome. 
Additionally, the resulting Na* current behaves contrarily to expectations in response 
to the astrocyte’s influence. Future work on our model should address our concerns 
regarding this particular parameter choice and the numerical stability of the model. 
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A 
Chapter 6 updaies | 
Inducing Flow Instabilities in Aneurysm 
Geometries via the Reynolds-Orr Method 


Alessandro Contri, Christina Taylor, Justin Tso, and Ingeborg Gjerde 


Abstract Ruptured intracranial aneurysms are the leading cause of hemorrhagic 
strokes. Although intracranial aneurysms are prevalent in as much as 4% of the 
adult population, most aneurysms do not rupture, and their growth and risk of 
rupture remains difficult to predict. Previous studies have identified abnormal wall 
shear stress patterns and blood pressure spikes, both features of unstable flows, as 
key factors in aneurysm behavior. While computational fluid dynamics has been 
enlisted to help study risk of rupture, no consensus has been established on how to 
impose unstable flow conditions. Here, we use Reynolds-Orr instability analysis to 
calculate flow perturbations that are capable of inducing unstable flow in 3D arterial 
geometries with and without aneurysms. These perturbations lay the foundation for 
future wall shear stress studies by providing mathematically consistent conditions 
for time-dependent Navier-Stokes simulations. 


6.1 Introduction 


Intracranial aneurysms are relatively common in adults, with a global prevalence of 
2-3%. Intracranial aneurysm rupture leads to aneurysmal subarachnoid hemorrhage, 
which is fatal in about 60% of cases and carries major risk of brain damage and 
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disability for survivors [1]. However, only 1-2% of intracranial aneurysms rupture 
annually, and most are asymptomatic [2]. The decision to treat unruptured aneurysms 
must therefore be carefully weighed against the risk of rupture, a metric that remains 
difficult to quantify. 

Aneurysm growth and rupture have been linked to local hemodynamic conditions, 
with wall shear stress (WSS) primary among them [3]. Computational fluid dynamics 
(CFD) studies allow for the high-resolution investigation of WSS in patient-specific 
arterial geometries through simulations of the Navier-Stokes equations. Several CFD 
studies have investigated the role of WSS in aneurysm rupture with conflicting 
results. While some studies have attributed risk of rupture to high WSS at thin and 
hard aneurysm walls, more recent studies have found strong relationships with low 
WSS at aneurysm rupture points (reviewed in [4]). Meng et al. have noted that both 
high and low WSS may increase rupture risk: high WSS causes degeneration of 
the cellular matrix and cell apoptosis, while low WSS and recirculation promotes 
an inflammatory environment that weakens the arterial wall and drives aneurysm 
growth [5]. 

Several studies have found rapid fluctuations in blood pressure to be a common 
risk factor for aneurysm rupture. Wouter ef al. found that 43% of ruptures occurred 
during stressful events, while Vlak et al. found high population-attributable risks 
for coffee drinking and vigorous physical exercise [6, 7]. Matsuda et al. also found 
high incidence of rupture in activities associated with the Valsalva maneuver, which 
results in sudden pressure changes across the aneurysm wall [8]. These perturbations 
in flow are highly relevant to aneurysm rupture, but are not accounted for in most 
CFD studies. 

These flow instabilities can be investigated using definitions of instability based 
on kinetic energy, which originated with Reynolds and Orr [9]. Analysis using these 
definitions, as revisited by Scott, provides an exact relation between a base flow 
and the evolution of the kinetic energy of a perturbation. The method is nonlinear 
with no approximations and yields a well-posed linear symmetric eigenvalue prob- 
lem that can be solved with standard methods. In the resulting eigenpairs the flow 
perturbations are captured in the eigenvectors, while negative eigenvalues indicate 
the potential for instability with respect to kinetic energy. In unstable eigenpairs, 
the magnitude of the eigenvalues indicates the initial growth rates of the perturbed 
kinetic energy. However, the most energetically unstable eigenpair is a function of 
both the magnitude of its eigenvalue and its eigenvector [9]. 

While the Reynolds-Orr eigenpairs may not exactly reflect physiological flow 
instabilities, they provide a reproducible way to induce flow instability. The readily 
available mathematical framework behind this approach allows this analysis to easily 
be reproduced in different arterial geometries. In this work, we investigate some of 
the nuances of using this approach to induce unstable flow on arterial geometries 
with and without aneurysms taken from the Aneurisk data set from Emory University 
[10]. 

In the remainder of this paper, we will first discuss the Reynolds-Orr method in 
more detail in Section 6.2.1 followed by the details of our time-dependent Navier- 
Stokes solver in Section 6.2.2. Next, we will present some of our initial findings on 
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using the Reynolds-Orr to induce instability in Section 6.3 before concluding with a 
discussion on our findings in Section 6.4. 


6.2 Methods 


In this section we consider the mathematical equations and the subsequent imple- 
mentation of the problem at hand, part of which was previously presented in [11]. 
The artery is considered as a pipe with a no-slip wall and prescribed pressures at 
the inlet and outlets. The absolute pressure is disregarded as it does not change the 
relative pressure drop in the model. 


6.2.1 Reynolds-Orr Instability 


The derivation of the eigenvalue problem from which we hope to detect the most 
unstable modes, given the base flow, follows the one presented in [9]. The derivation 
is nontrivial for pressure boundary conditions; therefore, we consider here only 
Dirichlet boundary conditions for the velocity. Let (u,p) be the solution to the 
following unsteady Navier-Stokes equations: 


0u-—vAut+u-Vu+Vp=f in Qx(0,T], (6.1) 
V-u=0 in Qx(0,T], (6.2) 

u=g on 0Qx(0,T], (6.3) 

u=u’ on Qx{0}. (6.4) 


Let now the couple (w,q) solve the same equations (6.4) with the same boundary 
condition g, but different initial data w°, such that w? + u°. Define (v,o) as the 
difference between the two solutions: 


v=u-W, O=p-gq. (6.5) 


Then, taking the two versions of (6.4) with the different initial data and subtracting 
them, (v, 0) satisfy: 


0,v-vAv+(u-Vu-w:-Vw)+Vo=90 in Qx(0,T], (6.6) 
V-v=0 in Qx(0,T], (6.7) 

v=0 on 09x(0,7], (6.8) 

v=v? =u°-w® on Qx {0}. (6.9) 


Following [9], the nonlinear term is rewritten as 
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u-Vu-—w:Vw=u-:Vu-u: Vw+u:Vww: Vw (6.10) 
=u-Vv+v:Vw=u-:Vv+v:Vu-v-Vv (6.11) 


Multiplying (6.9) by v, using (6.11) and integrating over Q yields 
(0,v,v) +a(v, v) +c(U, V, Vv) +c(v, U, Vv) —c(v, V, Vv) — b(v, 0) =0. (6.12) 


where (-,-) denotes the standard L? (Q) inner product and 


a(u,v) = v f vai, (6.13) 
Q 

c(u,v,w) = fvw, (6.14) 
Q 

b(u, p) = f v-up. (6.15) 
Q 


Defining now the following space: 
V:={ve H'(Q)?:V-v=0 and v=0 on dQ} (6.16) 
it can be observed from [9] that for any w € V one have 
c(w,v,v)=0, b(v,o) =0. (6.17) 


(6.12) can then be rewritten in the form 


1ć 2 2 
= y| = -y Vvi- - v-V "V .1 


1 
= -v f we -5 f Y 0u+vu. (6.19) 
Q 2 Jo 


The first term of the equation can be seen as the time derivative of the kinetic energy 
of the flow and the right hand side describes its evolution. Following standard 
arguments (i.e. integration in time of the left hand side etc.) we can say that the flow 
u is energy unstable at time ¢ = 0 if there exists a vo € V such that 


1 
-5 [vb (Gu0+ vvo -v | Vvo >o. (6.20) 
2 Jo Q 


At this point we define 


_ Bulv,v) 
Y a(v,v) 


1 
with Bu(v,w) = 5 J v! (Vu+ Vu! )w (6.21) 
Q 


Comparing (6.19) with (6.21) we can reformulate an equivalent instability condition. 
Namely, the flow is unstable at t = 0 if there is a vo € V such that 


Av <1. (6.22) 
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Given that a(v,v) = v||Vv||?, we have in particular that if Sv such that 4 = 
infozyey Ay < —1 and thus 


1 


ð 2 
(rim = 2v(1+Ay) > 0, (6.23) 


the flow is unstable. Recalling Poincaré’s inequality, stating ||v||7 < Cp||Vv||?7,_ Yve 
V, we can deduce 


Ot 


We conclude that the most negative A, leads to the most unstable mode. However, as 
noted in [9], the most unstable mode may not be the most persistent. Additionally, 
[9] proves that the solution to A = infozyey Ay solves the eigenvalue problem (A, v) € 
(R, V) such that 


we f we) > SN AEAN (6.24) 
Q Cp 


Bu(v,w) =Aa(v,w) YweV. (6.25) 


The detailed derivation can be found in [9]. 


6.2.2 Unsteady Navier-Stokes FEM Discretization 


We looked to study the evolution in time of the potentially unsteady modes detected 
using the method described previously. In order to progress in time we needed to 
implement a time-dependent Navier-Stokes solver. To do so, we used a splitting 
method, which despite its low accuracy in time is fast and easy to implement from 
the computational point of view. The motivation behind splitting schemes is to 
consider the first two equations in 6.4 separately. We chose a modified version of 
Chorin’s method [12], the so-called incremental pressure correction scheme (IPCS) 
[13], which gives improved accuracy compared to the original scheme at little extra 
cost. 

The IPCS scheme involves three steps. First, we compute a tentative velocity u* 
by advancing the momentum equation (6.4) by a midpoint finite difference scheme 
in time, using the pressure p” from the previous time interval. We also linearize the 
nonlinear convective term by using the known velocity u” from the previous time 
step: u” - Vu". The variational problem for this first step is 


((u* —u”)At,v) + (u” - Vu", v) +.a(u"*?, v) — b(v, p”) (6.26) 


+(p"n, Vag — v(Vu"*? -n, v)ao = (£"*!,v) (6.27) 


where 


vas f uvds, ow? = (u” +u”*!)/2. (6.28) 
dQ 
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Since we are modeling flow in a ’pipe” (the artery) with known inflow and outflow, 
the scheme has the hidden assumption that the derivative of the velocity in the 
direction of the channel is zero at the inflow and outflow, corresponding to a flow 
that is ”fully developed,” or that the flow field doesn’t change significantly upstream 
or downstream of the domain. 
The second step is to use the computed tentative velocity to compute the new 
pressure p”: 
(Vp""!,Vq) =(Vp",Vq)- At! (V -u*,q) (6.29) 


Finally, we compute the corrected velocity u”+!:; 


(ut! vy) = (u*, v) — Ar(V(p""! — p”),v). (6.30) 
In summary, we solve the incompressible Navier-Stokes equations efficiently by 


solving a sequence of three linear variational problems in each time step. The FEM 


space chosen to solve the problem is (uy, pp) € (V2. m v!) where, if we let vi denote 


the space of C* piecewise polynomials of degree k on a regular mesh of the domain 
Q, the vector valued polynomial space is 


VĚ „= {v € (Vi)? : Vlaa\r = 0}. (6.31) 


6.3 Results 


For our numerical experiments we used a single artery mesh both with and without 
a aneurysm. These geometries are shown in Figure 6.1. 


Case 0 Case 1 


Fig. 6.1: The artery geometry without (Case 0) and with (Case 1) an aneurysm used 
for our numerical experiments. 
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As previously stated, eigenvalues with negative sign and magnitude greater than 
one are necessary to induce energy instability. Eigenpairs were computed using a 
shifted power iteration with LU preconditioning to target eigenvalues near -1. We 
observed that without a sufficient pressure drop, eigenvalues less than or equal to -1 
were not guaranteed. Using Ap = 600Pa (4.5mmHg) and 1000Pa (7.5mmHg) both 
yielded at least two eigenvalues less than -1 in both the healthy artery and terminal 
aneurysm. Tables 6.1 and 6.2 show the eigenvalues reported by the solver. Figures 
6.2 and 6.3 show example perturbations associated with the Ap = 600Pa pressure 
drop case. 


Table 6.1: Eigenvalues for the healthy artery without an aneurysm for different 
pressure drop values. 


N Ap=600Pa Ap=1000Pa 
(4.5mmHg) (7.5mmHg) 

0 -1.116 -1.778 

1 -1.091 -1.623 

2 -0.874 -1.319 

3 -0.775 -1.163 

4 0.871 1.265 

5 0.964 1.402 


Table 6.2: Eigenvalues for the artery with an aneurysm for different pressure drop 
values. 


N Ap=600Pa Ap=1000Pa 
(4.5mmHg) (7.5mmHg) 

0 -1.154 -1.865 

1 -1.125 -1.697 

2 -0.903 -1.404 

3 -0.785 -1.191 

4 0.875 1.160 

5 0.970 1.279 

6 s 1.443 


* The 600Pa pressure drop did not yield a sixth eigenvalue near +1. 


With the eigenpairs in hand, it remained to simulate the kinetic energy evolution 
of the eigenpairs with eigenvalues less than -1. The Navier-Stokes equations are 
notoriously difficult to simulate especially in the face of physical instabilities. In the 
case of our scheme a severely prohibitive time step was needed to ensure numerical 
stability. The prohibitive time step coupled with a lack of computing resources 
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Fig. 6.2: Example of the base flow without and with perturbation added in the healthy 
artery case for Ap = 600Pa. 
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Fig. 6.3: Example of base flow without and with perturbation added in the aneurysm 
case for Ap = 600Pa. 


restricted us from being able to perform in-depth analysis of the kinetic energy 
evolution. However, we were able to calculate the kinetic energy evolution for the 
above eigenpairs over a very short time period. Figures 6.4 and 6.5 show the L2 
norm of the velocity (the kinetic energy) of the solutions over a time period of 0.1 
seconds. 


6.4 Discussion 


From our initial results, only four eigenpairs seem to induce unstable flow as indicated 
by increasing kinetic energy despite there being 16 negative eigenvalues. However, 
as shown in Figure 6.4 for the 1000Pa case with A = —1.319, it seems possible 
for the kinetic energy to initially decrease before increasing. This is not completely 
unexpected as [9] showed the kinetic energy can oscillate, though in their experiments 
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Fig. 6.4: Evolution of the L2 norm of the velocity in time for healthy artery eigenpairs. 


this oscillation was coupled to an overall decay after an initial increase in energy. 
Our short simulation time may not have been sufficient to capture instabilities that 
require more time to develop. Additionally, as was noted in [9] a larger magnitude 
negative eigenvalue does not necessarily guarantee the most unstable mode due to 
the instability’s additional dependence on the gradient of the perturbed velocity. In 
Figure 6.5 this can be seen in the 1000Pa case by the eigenpair with 4 = —0.785 
inducing an initial increase in energy while the eigenpair with A = —1.154 did not. 
While the restrictive time step required by the Navier-Stokes solver did not allow 
us to complete all initial goals of this project, we have been able to show some 
preliminary results showing that perturbations computed using the Reynolds-Orr 
method can induce energy instability for a given geometry. Future work would 
entail longer simulations to confirm whether all negative eigenvalues with sufficient 
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Fig. 6.5: Evolution of the L2 norm of the velocity in time for the aneurysm eigenpairs. 


magnitude do lead to instability eventually. Additionally, we have completed the 
work started in [11] by correcting errors that had previously made when calculated 
the eigenmodes. The codes produced here provide the framework needed to continue 
studying aneurysm rupture risk through the Reynolds-Orr instability analysis when 
provided sufficient computing resources. 


6 Inducing Flow Instabilities in Aneurysms 89 


References 


10. 


12. 


13. 


Marilyn M Rymer. Hemorrhagic stroke: intracerebral hemorrhage. Missouri medicine, 
108(1):50, 2011. 

Gabriel JE Rinkel, Mamuka Djibuti, Ale Algra, and J Van Gijn. Prevalence and risk of rupture 
of intracranial aneurysms: a systematic review. Stroke, 29(1):251-256, 1998. 

Mahsa Dabagh, Priya Nair, John Gounley, David Frakes, L Fernando Gonzalez, and Amanda 
Randles. Hemodynamic and morphological characteristics of a growing cerebral aneurysm. 
Neurosurgical focus, 47(1):E13, 2019. 

Yuichi Murayama, Soichiro Fujimura, Tomoaki Suzuki, and Hiroyuki Takao. Computational 
fluid dynamics as a risk assessment tool for aneurysm rupture. Neurosurgical focus, 47(1):E12, 
2019. 

H Meng, VM Tutino, J Xiang, and A Siddiqui. High wss or low wss? complex interactions of 
hemodynamics with intracranial aneurysm initiation, growth, and rupture: toward a unifying 
hypothesis. American Journal of Neuroradiology, 35(7):1254—1262, 2014. 

Wouter I Schievink, John M Karemaker, Loes M Hageman, and Dries JM van der Werf. Circum- 
stances surrounding aneurysmal subarachnoid hemorrhage. Surgical neurology, 32(4):266- 
272, 1989. 

Monique HM Vlak, Gabriel JE Rinkel, Paut Greebe, Johanna G van der Bom, and Ale Algra. 
Trigger factors and their attributable risk for rupture of intracranial aneurysms: a case-crossover 
study. Stroke, 42(7):1878-1882, 2011. 

Masayuki Matsuda, Kazuyoshi Watanabe, Akira Saito, Ken-ichi Matsumura, and Masaharu 
Ichikawa. Circumstances, activities, and events precipitating aneurysmal subarachnoid hem- 
orrhage. Journal of Stroke and Cerebrovascular Diseases, 16(1):25—29, 2007. 

L Ridgway Scott. Kinetic energy flow instability with application to couette flow. Technical 
report, Technical Report TR-2020-07, 2020. 

AW Bergersen, HA Kjeldsberg, and K Valen-Sendstad. Kvslab. http://ecm2.mathcs. 
emory .edu/aneuriskweb/index, 2021. Accessed August 2022. 


. Nanna Berre, Gabriela Castro, Henrik Kjeldsberg, Rami Masri, and Ingeborg Gjerde. A 


computational study of flow instabilities in aneurysms. In Simula SpringerBriefs on Computing, 
pages 63-75. Springer, Cham, 2022. 

Alexandre Joel Chorin. The numerical solution of the navier-stokes equations for an incom- 
pressible fluid. Bulletin of the American Mathematical Society, 73(6):928—931, 1967. 
Katuhiko Goda. A multistep technique with implicit difference schemes for calculating two-or 
three-dimensional cavity flows. Journal of computational physics, 30(1):76—95, 1979. 


Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International 
License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and 
reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the 
source, provide a link to the Creative Commons license and indicate if changes were made. 


The images or other third party material in this chapter are included in the chapter’s Creative Commons 


license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter’s 
Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the 
permitted use, you will need to obtain permission directly from the copyright holder. 


ÁA 
Chapter 7 ‘pdates 
Impact of Pathological Vascular Remodelling on 
Right Ventricular Mechanics 
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Abstract Pulmonary arterial hypertension (PAH) is a rare disorder characterized by 
elevated blood pressure and pulmonary vascular resistance, often followed by right 
ventricular hypertrophy and heart failure. The effect of PAH and its treatments on the 
mechanics, function, and remodelling of the right ventricle (RV) is currently not well 
understood. To study cardiac biomechanics and functionality as PAH progresses, 
we implemented a computational model of the heart simulating right ventricular 
maladaptive remodelling. Our Windkessel-based model, which accounts for direct 
ventricular interaction and the presence of the pericardium, is utilized to simulate 
various disease stages of PAH. We find that the pericardium has a larger effect on 
heart performance than ventricular interaction through the septum. We also examined 
the effectiveness of two treatments, ventricular assist device (RVAD) and atrial 
septostomy, on diseased hearts. We show that while both pulsatile and continuous 
RVADs restore cardiac function, pulsatile RVAD improves cardiac output 29.4% 
more than continuous RVAD. We also demonstrate that atrial septostomy improves 
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cardiac output by 19.5%. Our model can be further extended by simulating the heart’s 
response to other treatments such as extracorporeal membrane oxygenation (ECMO), 
and by incorporating ventricular remodelling growth simulations and finite-element 
ventricular modelling. 


7.1 Introduction 


Pulmonary arterial hypertension (PAH) is characterized by high blood pressure and 
high pulmonary vascular resistance in the arteries in the lungs [1]. As the pressure 
and resistance in the pulmonary arteries and capillaries rises, the right ventricular 
afterload is elevated. Over time, there is a larger demand for work done by the heart 
to achieve the same level of cardiac output as that of a healthy one. To increase 
the amount of contractile force generated by the heart, the right ventricle (RV) 
undergoes abnormal enlargement and increases its wall thickness, stiffness, and 
contractility. This phenomenon, which is observed in almost all PAH patients, is 
referred to as right ventricular hypertrophy [2]. Moreover, the heart compensates for 
the elevation in pulmonary arterial pressure and resistance via vascular proliferation 
and remodelling of small pulmonary arteries [3]. In advanced stages of PAH, the 
severe dilation of the RV causes septal bowing, i.e., a leftward shift of the septum 
that lowers left ventricular volume and limits left ventricular filling. Hence, left 
ventricular output is reduced and, in turn, cardiac output and mean arterial pressure 
decrease [4]. All these compensatory mechanisms enable the heart to maintain its 
performance in the early stages of PAH; however, sustained structural and functional 
changes in the RV can eventually become detrimental to cardiac function and may 
even result in right heart failure and death of a patient [5]. 

Although life-threatening, PAH remains incurable. Treatment strategies are avail- 
able to alleviate the vascular symptoms (e.g., vasoconstriction) [6], and to delay 
disease progression by unloading the RV. Therapies proven to be effective in PAH 
patients include: drugs in the form of vasodilators, inotropes and anticoagulants; 
extracorporeal membrane oxygenation (ECMO) devices; mechanical support via a 
right ventricular assist device (RVAD); atrial septostomy; or lung transplantation. 
[7]. Multiple in situ and in silico investigations focus on PAH, yet there is still a 
gap in our understanding of the biomechanic and hemodynamic effect of PAH on 
the RV and its subsequent remodeling (and failure). [3, 4, 8, 9, 10] Furthermore, 
the progressive remodeling of the RV is inevitably influenced by its surrounding 
environment; in other words, the septum, the pericardium, the left ventricle (LV), 
and baroreflex must be considered in order to investigate RV performance in PAH 
as accurately as possible. Some computational studies have already looked into 
changes in heart function with ventricular interdependence under the influence of 
pericardium and baroreflex. [11, 12, 13] Nonetheless, most rigorous biomechanical 
and hemodynamic models have yet to be utilized for research involving treatments 
of PAH and right ventricular remodelling. Our objective was to develop a realistic 
computational model that simulates the heart conditions in a PAH setting to gain 
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insight into the various maladaptations characteristic of this condition as well as the 
effect of PAH treatments on RV mechanics. Our theoretical analysis supplements 
experimental research techniques to further enhance our knowledge of PAH ther- 
apies and how the progressive adverse remodeling impacts the biomechanics and 
hemodynamics of the heart. 


7.2 Methods 


In this work, we implement a closed-loop, lumped-parameter model of the heart and 
the circulatory system [8, 12] to study the impact of increased pulmonary vascular 
resistance on RV hemodynamics. This Windkessel-based model captures the effects 
of PAH at different severities while accounting for direct ventricular interaction via 
the interventricular septum, the pericardium constraint, baroreflex, and the heart’s 
dynamic growth and remodelling in response to varying loads. Moreover, our model 
is utilized to examine two treatment options for PAH-induced RV dysfunction. In par- 
ticular, we investigate the impact of RVAD and atrial septostomy on cardiovascular 
hemodynamics in the PAH setting. 


7.2.1 Base Model of PAH 


Our model of the cardiovascular system combines elements of the cardiac anatomy 
with those of the systemic and pulmonary circulations. The pressure, volume, and 
blood flow through the circulatory system are represented using a set of ordinary 
differential equations. The relationship between the pressure and volume for each 
component in the circulatory system can be represented linearly by 


V=CP+V4q, 


where V is volume, P is pressure, C is the compliance of the vessel, and Vg is the 
volume of the vessel at zero pressure. While C and V4 of the arteries and veins are 
typically constant (unless there are physiological or pathological changes), C and 
Va of the heart chambers are time-varying periodic functions that represent active 
contraction [14]. The circulatory system model consists of individual elastic cham- 
bers and a series of resistors and diodes with fixed and variable capacitors branching 
from the main circuit (Figure 7.1) [15]. The systemic and pulmonary vasculatures 
are divided into their corresponding resistive and fixed capacitive components due to 
their elastic nature. Heart valves reinforce uni-directional blood flow, so they are rep- 
resented as diodes. In addition, the mitral and tricuspid valves exhibit non-negligible 
resistive elements to block blood flow between heart chambers. Systole and diastole, 
the two phases of the cardiac cycle, are captured by variable capacitors that represent 
the heart chambers. 
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Fig. 7.1: Base Windkessel circulation model describing the coupling between the 
heart and the systemic and pulmonary circulation, including RVAD and atrial sep- 
tostomy treatments. 


The volume of each individual compartment is related to the rate of inflow and 


outflow Vla 
ETE Qin sF Qout, ] -1 


where V; is the volume of chamber i, and Qin and Qout represent the flow in and 
out of the chamber, respectively. We approximate the blood flow at each component, 
which is driven by pressure difference, to be linearly dependent on pressure through 
the following relation: 


Piss -P : 

pstream downstream 

Q = —— (7.2) 

R 

where R is the resistance of the given connection. 
Ventricular and atrial pumping characteristics are modelled by modifying time- 

varying elastance, e(t). Mathematically, instantaneous ventricular pressure P(t) can 


be related to instantaneous ventricular volume V(t) via 
P(t) = [Pes(V) F Pea(V)] e(t) +Pea(V), (7.3) 


where 
Pes(V) = Ees(V — Vo), (7.4) 


Pea(V) = a(ef V- — 1), (7.5) 
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h|sin(r-)+1], ifo<r< Ye 
e(t) = (7.6) 
0.5exp | (« - ze) ir}, otherwise. 


Here, Pes is end-systolic pressure, Peq is end-diastolic pressure, Ees is end-systolic 
elastance, Tes is time required to achieve end-systole, r is the time constant of 
relaxation, and Vo, œ, 8, are constant model parameters in the end-diastolic pressure- 
volume relationship. The parameters of the base model were adapted, based on 
Punnoose et al. [8], to capture different stages of PAH (mild, moderate, severe, 
and very severe, i.e., cardiogenic shock). We initialized the model with standard 
parameters for a healthy heart, and then varied its parameters to match those of 
the corresponding phase. When the parameters were altered for each diseased state, 
common RV function metrics, such as the stroke volume, the ejection fraction and 
cardiac output, changed too. We could thus observe a distinct pressure-volume 
relationship for each degree of right-sided heart failure before implementing any 
treatments. Overall, the most relevant parameters in the circulation model include 
systemic and pulmonary vessel resistance, compliance and elastance. The varying 
effects of most of these parameters on the model’s behaviour was further explored 
by performing a sensitivity analysis. 


7.2.2 Ventricular Interaction and the Pericardium 


In this study, we aimed to improve the realism of the previously described lumped- 
parameter model by accounting for the interventricular interaction via the septum 
wall. The left and the right ventricles act independently of each other; their time- 
varying volume-elastance relation is represented by the system shown above in 
Equation 7.6. Following Santamore et al. [11], three time-varying elastances (right 
ventricular free wall E,„,p, left ventricular free wall E, f, and septum wall Es) were 
considered, so that the interdependence between the two ventricles was simulated 
by our model. The pressures of both ventricles at the end of diastole (Peq) becomes 


Peary (V) = ale" a1), (7.7) 
where 
By = Biot+mriVry, (7.8) 
and 
Pearv(V) = a(e?” VW) _ 1), (7.9) 


where 
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By = Bro + mir Viv (7.10) 


Bio and B,o represent the ventricular coefficients when the volume is equal to zero. 
my, is the sensitivity of Aj, to changes in right ventricular volume while m1, is the 
sensitivity of A,,, to changes in left ventricular volume. The pressure development in 
each ventricle during systole is a function of pressure in the opposite ventricle. For 
the end-systolic relation, the pressure at the end of systole (Pes) of both ventricles is 
formulated as 


Es Ewf Ew ¢ Pry 
P = — — (Vp — Vo) + ————, 7.11 
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Fig. 7.2: Circuit used to simulate cardiovascular system with interaction from San- 
tamore et al.[11] 


We then proceeded to extend our model to account for the presence of the peri- 
cardium, a fibrous sac that encloses the four chambers of the heart. We followed the 
implementation from Burkoff et al. [12] in which the pericardial pressure (P pericard) 
is assumed to be equal to an exponential function of the sum of instantaneous LV 
and RV volumes: P pericard = @p exp[Bp (Viv +V,v)]. 
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7.2.3 Right Ventricular Remodelling and Baroreflex 


The heart continuously regulates its output to match the needs of the body, but it 
also adapts its anatomy and behavior to changes in pulmonary arterial pressure. 
During PAH, the heart undergoes long-term deformations in the form of an increase 
in contractility, and subsequent augment in right ventricular wall volume, Vai. 
Altered contractility, which is associated with changes in sarcomere length from the 
ideal optimal value, is evaluated in our model using the formulation from Arts et al. 
[16]. The sarcomere length can be calculated as 


1/3 
A= (13 =] , (7.13) 


wall 


and the change in contractility required to reach such sarcomere length [17] is 


computed as 
l+a 


= Dga ebo Loma)’ (1.14) 


where a = 0.2, b= 4um=! [16], Lo is the desired sarcomere length, and Ls max is the 
maximum sarcomere length throughout the cycle. Once C has been calculated, wall 
volume (Vai) and right ventricular contractility (Ees,rv) are respectively updated 
in an iterative process. Mathematically, a linear approximation from the previous 
iteration is used in both cases: Egy „y = Cea, and V) = cl. 

A more short-term cardiac adaptation to changes in heart rate and contractility is 
baroreflex control. This homeostatic mechanism increases the heart rate when arterial 
pressure decreases, and vice versa [18]. Changes in the heart rate were simulated 
using the linear approximation [13] T = Gparo(Pao — 120) + 0.855, where T, is the 
cardiac period, Pao is the systolic aortic pressure and G paro (the baroreflex gain) is 
fixed at 0.005 s/mmHg. 


7.2.4 Atrial Septostomy 


Atrial septostomy is a palliative treatment for medically refractory PAH. The pro- 
cedure involves using a catheter to intentionally create a hole, called an atrial sep- 
tostomy defect (ASD), in the interatrial septum. The direct blood flow from the right 
atrium into the left atrium increases cardiac output, thus alleviating right atrium 
pressure. As a result, the overall outflow resistance and afterload experienced by 
the right-side of the heart is reduced, enhancing survival. We account for atrial 
septostomy in our model by introducing a new lump-parameter for septostomy re- 
sistance, Rs¢ pr, and adding blood flow through the interatrial septal orifice governed 
by Qsept = K Ay|AP|, where A = zR? is the area of the septostomy in cm?, AP is 
the pressure gradient across the orifice in mmHg, and K = 2.66. 
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7.2.5 Right Ventricular Assistive Device (RVAD) 


PAH often leads to right ventricular failure, so often assistive devices (RVADs) are 
required. These devices can be used in a continuous flow capacity, where the flow 
continuously provides support to the RV independent of time; or as a pulsating 
RVAD, where the flow is time dependent and in sync with the cardiac cycle. Both 
continuous and pulsatile flow pumps have previously been tested computationally 
and experimentally [9, 19]. Additionally, the RVAD can be implemented in either 
the right ventricle or right atrium. To implement the RVAD into the model, we add 
an additional flow, Q,„yvaa. The flow into the pulmonary artery will then be given by 


dV pa 
dt = Opna-QpvtQrvad- (7.15) 


RVAD will take fluid from the right ventricle or atrium, which results in either 


dV -y 

dt = Orr — Opwy — Orvad (7.16) 
or dy 

= sv ` Vtr — Urvad> 7.17 

TTE = Osv — Qir — Orvad (7.17) 


respectively. The form of Q;,aq depends on the type of RVAD. The continuous 
flow RVAD only depends on the difference in pressure between the source (the right 
atrium or ventricle) and the pulmonary artery. Analytically, this can be approximated 
using a linear relationship [8], where Q;-).aq is time independent. We also implement 
the pulsatile flow pump as described in Gohean et al. [9], which has demonstrated 
better performance compared to the continuous flow pumps [19], when properly 
synchronized. Now Q+vaa is time dependent and follows the formulation in Gohean 
et al. [9], given as 


2V, 2r (t-Ta i 
T (3-4 cos (2=6772)) if t € [Ta,Ta +T], 


Qrvaa(t, AP) = f (7.18) 


otherwise. 


where T4 is the time delay from the start of the cardiac cycle, T; is the duration of 
the pulse, and V, is the stroke volume in mL. The time, t, is normalized to the start 
of the cardiac cycle. 


7.3 Results 


We successfully implemented a ventricular hemodynamic model and simulated PAH 
based on data from Punnoose et al. [8] to study the effects of this aggressive dis- 
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ease on RV hemodynamics. The right and left ventricular hemodynamic parameters 
reproduced by our model are reported in Table 7.1. 


Table 7.1: Stroke volume (SV), cardiac output (CO), and ejection fraction (EF) from 
the model for healthy, severe PAH, continuous RVAD treatment, and pulsatile RVAD 
treatment. Stroke volume is given in ml and cardiac output in L/min. 


Left ventricle Right ventricle 
Condition SV CO EF SV CO EF 
Healthy 72.722 4.363 0.717 72.722 4.363 0.687 
Severe PAH 33.003 1.980 0.617 33.002 1.980 0.209 
Continuous RVAD 46.594 2.795 0.646 22.840 1.370 0.155 
Pulsatile RVAD 56.302 3.378 0.658 19.901 1.194 0.152 


An important indicator of the severity of PAH is RV function and morphology 
[20]. Figure 7.3 shows pressure-volume loops (PV loops) corresponding to various 
severities of PAH. The shapes of the PV loops, which provide information on cardiac 
performance, load and coupling, change with progression of PAH. To simulate the 
different stages of PAH, the elastance, resistance and compliance parameters as well 
as the heart rate were modified based on values from previous studies [8]. As the 
disease progresses, resistance in the pulmonary artery increases, heart rate increases, 
and end-systolic elastance increases in the left ventricle, while it decreases in the 
right ventricle. For the LV, the normal shape of the PV loop is rectangular, whereas 
it is usually rounded for the RV. For both ventricles, the PV loops become more 
narrow, which is indicative of reduced stroke volume. 
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Fig. 7.3: PV loops for healthy and diseased hearts, from mild PAH to PAH-induced 
heart failure. Dark blue represents the most severe case, cardiogenic shock (CGS). 
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7.3.1 Sensitivity Analysis 


A sensitivity analysis was a useful tool to determine the correlation between input 
model parameters and the resulting hemodynamic features of the system. By altering 
each model parameter, one at atime, we were able to better interpret the model outputs 
when the inputs were changed in PAH treatment simulations. The model behavior 
changed drastically when the pulmonary arterial resistance (Rpa) was altered, as 
illustrated in Figure 7.4. The Rpa values plotted correspond to 0.25, 0.50, 0.75 
and 1 mmHg.s/mL. As Rpa increased (as in PAH), the stroke volume decreased 
substantially. The total mechanical energy generated by ventricular contraction also 
decreased, since the area encompassed within the PV loops (i.e., the stroke work) 
decreased with higher R pa values. In addition to pulmonary arterial resistance, other 
parameters that also had a considerably high impact on the model outputs were RV 
end-systolic elastance and the exponential parameter (£) for the right ventricular 
end-diastolic PVR (7.5). 
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Fig. 7.4: PV loops for several degrees of pulmonary arterial resistance. 


7.3.2 Ventricular Interaction and the Pericardium 


As explained in subsection 7.2.2, we aimed to improve the realism of the model by 
including the interaction via the interventricular septum. We implemented three time- 
varying elastances with the same values used by Santamore et al. [11], and simulated 
different degrees of PAH by changing the values of resistance as in Punnoose et al. 
[8], and keeping the same values of elastance from Santamore et al. [11]. Figure 7.5 
shows the different pressure-volume (PV) loops of the left and right ventricles when 
the interventricular interaction is turned on. The overall behavior of the PV loops is 
the same as those without interaction: there is a leftward shift of the LV PV loop anda 
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rightward shift of the RV PV loop with increasing PAH severity. This is likely due to 
the RV expanding and pushing against the LV, causing the right ventricular volume to 
increase and the left ventricular space to decrease, eventually leading to ventricular 
septum bowing. Compared to Figure 7.3, Figure 7.5 displays PV loops that shift 
more uniformly along an axis as severity increases. One possible explanation for 
this behaviour is the use of the same elastance values for the ventricular interaction 
simulations at each severity level. This behaviour, however, is not as strongly evident 
as in the previous simulations, which could be explained by the fact that we kept 
the same elastance from Santamore et al. [11] for each PAH degree, and because 
we used a combination of two parameter sets (Punnoose ef al.[8] and Santamore 
[11]). We can still conclude that the introduction of the interventricular interaction 
does not lead to strong changes for this model. This result is in agreement with other 
studies that report small effects from the septum [11]. 
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Fig. 7.5: Pressure-Volume loops with different degrees of PAH when the ventricular 
interaction via the septum is included. Dark blue represents cardiogenic shock (CGS). 


Secondly, the presence of the pericardium represents an additional constraint to 
both ventricles. Figure 7.6 shows the effect of the pericardium in a healthy PV loop. 
When this element was added to our model, it lead to a reduction in the overall left 
ventricular stroke volume and an increase in the end-diastolic pressure of the right 
ventricle. The PV loop shifts slightly upward off the end-diastolic pressure-volume 
relationship (EDPVR), which is representative of the heart’s ventricular passive 
compliance curve. In addition, the stroke volume is smaller and the peak pressure 
is lower in presence of pericardium. The stroke volume reduction observed in our 
model is consistent with physiological and clinical observations that the pericardium 
constricts heart dilation in the diastolic phase. 
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Fig. 7.6: Pressure-Volume loops of healthy condition with and without pericardium 


7.3.3 Remodelling the Right Ventricle 


Including the remodelling effects of the RV resulted in improved performance of 
the heart. Computing the change in contractility, as described by Equation 7.14, 
yielded improved model output for an increase in pulmonary pressure (Figure 7.7). 
Additionally, the arterial baroreceptor reflex system changed the model behavior 
in presence of PAH. The heart rate and contractility (heart rate effect shown in 
Figure 7.8), which are two parameters that have a significant effect on the model’s 
output, are adjusted by the baroreflex to control fluctuations in aortic blood pressure. 
Thus, impaired baroreflex control of the heart in PAH would typically contribute to 
worse outcome of the disease. 
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Fig. 7.7: Effect of remodelling the RV on the pressure volume relationship in the left 
and right ventricle. 
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Fig. 7.8: Heart rate control, characteristic of baroreflex, has a high impact on the 
shapes of the PV loops. 


7.3.4 Atrial Septostomy 


Consistent with the left- and right-ventricular hemodynamics behavior observed in 
Punnoose et al., our atrial septostomy defect simulation with resulting right-to-left 
shunt yields minimal leftward shifts in both the right atrial and right ventricular 
pressure-volume (PV) loops toward lower volumes and pressures [8]. Simultane- 
ously, a rightward shift in the left atrial and left ventricular PV loops reflect increased 
filling pressures and end-diastolic volumes as shown in Figure 7.9. 

Based on the results from our model, atrial septostomy significantly increased 
left ventricular stroke volume as a result of elevated end-systolic pressure and ex- 
panded end-diastolic volume. Meanwhile, the left ventricular end-diastolic pressure 
and end-systolic volume experienced no significant change with atrial septostomy. 
Collectively, these results indicate atrial septostomy causes an increase in left ven- 
tricular cardiac output with negligible change in the systemic vascular resistance. 
The right ventricular pressure, volume, and stroke volume remained largely the same 
with or without atrial septostomy and are thus determined to be minimally responsive 
to atrial septostomy treatment. 

Prior to atrial septostomy treatment, the left and right atria are isolated, pressurized 
elastic chambers, with the right atrium having a higher pressure than the left atrium. 
The chamber pressure values are within reasonable ranges of 1 to 7 mmHg for the 
right atrium and 1 to 6 mmHg for the left atrium, as the atria are low pressure 
chambers meant for gathering blood from veins. After atrial septostomy treatment, 
however, we observed an increase of the end-systolic pressure and the end-diastolic 
volume in the left atrium as well as a slight decrease in the right atrium end-systolic 
pressure and end-diastolic volume. This is likely due to the two atria being connected 
via the atrial septal defect post-surgery, driving the pressures in the left and right 
atria to an equilibrium. 
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Our model observed a decrease in interatrial pressure gradient. This is demon- 
strated in Table 7.2 by a 34% elevation in left atrial end-diastolic pressure to match 
an 11% decrease in right atrial end-diastolic pressure. This outcome aligns with our 
prediction that the two atrial chambers will reach a pressure equilibrium when two 
heart chambers are connected via the atrial septal defect. 
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Fig. 7.9: Effect of atrial septostomy (ASD) on the pressure-volume relationship in 
each heart chamber. 


Table 7.2: Heart chamber end-diastolic pressures (EDP) with and without atrial 
septostomy. 


Heart Chamber EDP, no septostomy EDP, septostomy Pressure Ratio 
(mmHg) (mmHg) 

Right Atrium 3.142 2.795 0.8897 

Right Ventricle 2.237 2.046 0.9147 

Left Atrium 2.084 2.795 1.3415 

Left Ventricle 1.160 1.687 1.4550 


7.3.5 Right Ventricular Assistive Device 


We find that our RVAD implementation behaves similarly to the one in Punnoose et 
al. [8] and Gohean et al. [9] (Table 7.1 and Figure 7.10). Simulated PAH reduced left 
ventricular cardiac output, stroke volume and ejection fraction (Table 7.1). This is 
evident in Figure 7.10, with a leftward shift in the left ventricular PV loops indicating 
left ventricular unloading. This is improved by pulsatile and continuous RVAD at a 
flow rate of 3.5 L/min from the right atrium or ventricle, indicated by a rightward 
shift in the left ventricular PV loops for both treatments, with the pulsatile flow 
performing better than the continuous flow RVAD (Figure 7.10 right). Conversely, 
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the right ventricular end-systolic and diastolic pressures increase in the simulated 
PAH. As shown in Figure 7.10 (left), a rightward shift in the PAH PV loop is seen, 
which is not improved by either the pulsatile or continuous RVAD. 
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Fig. 7.10: Ventricular PV loops in healthy, severe PAH and treatment cases. The 
right ventricular PV loops for severe PAH shift to the left, but are partially recovered 
by pulsatile or continuous RVAD (left). Left ventricular PV loops in healthy, PAH 
and treated simulations of hemodynamics (right). 


7.4 Discussion 


In this study, we have modeled the influence of PAH on the mechanics of the 
circulatory system at various stages of the disease, based on the results from Punnoose 
etal. [8]. A sensitivity analysis demonstrated the significant dependence of our model 
on the resistance in the pulmonary artery, but many factors beyond resistance, such 
as contractility and stiffness, have an effect on the model too. To improve upon the 
base model, the addition of the interaction between the ventricular walls through the 
septum [11] and the effects of the pericardium [13] were included. We find a larger 
influence from the pericardium, which agrees with previous results [11, 13]. Further, 
we investigated long and short term adaptions of the heart to the increase in pressure 
of the pulmonary artery, as well as the effect of treatments to PAH. Specifically, the 
target of our computational study was to provide insight into how atrial septostomy 
and RVAD affect heart function [8, 9]. 

Our model exhibited hemodynamic effects of atrial septostomy similar to those 
observed in Punnoose ef al. for all four heart chambers. This validated our im- 
plementation to study the effect of atrial septostomy on PAH-induced right-side 
heart failure [8]. Results from our simulation indicate that atrial septostomy indeed 
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restores the heart to a healthier functional state under the assumption of ceteris 
paribus, consistent with clinical observations. The end-diastolic pressure-volume 
relationship (EDPVR) on the PV loop represents the passive compliance of the heart 
ventricles. With atrial septostomy, the EDPVR remained the same in both the left 
and right ventricles, implying that the heart ventricular tissue stiffness is not affected 
by the treatment. The end-systolic pressure-volume relationship (ESPVR), which 
represents ventricular tissue contractility, remained the same for the left ventricle 
and decreased for the right ventricle with atrial septostomy treatment, denoting less 
contractility. The introduction of atrial septostomy in the interatrial septum likely 
decreased the right ventricle’s need to contract to maintain the same cardiac output 
under PAH. Collectively, these data suggest that, while atrial septostomy does not 
change the myocardium intrinsic biomechanical properties, the surgical procedure 
seems to promote cardiac output, improve heart performance, and reduce heart fail- 
ure symptoms. Summed together, our model provided results consistent with current 
clinical understanding of atrial septostomy effects on cardiovascular biomechanics, 
physiology, and literature. 

The RVAD resulted in improved cardiac performance for a severe PAH setting 
by restoring cardiac output and ejection fraction in the LV. In our simulations, we 
observe better performance from the pulsatile flow compared to the continuous flow 
pumps. Despite the improved performance, the pulsatile flow pumps typically suffer 
from durability issues [9], so continuous pumps must be investigated too. Using this 
model, the influence of both types of pumps on the hemodynamics can continue to 
be studied in future research. 

Our model was able to capture the effects of PAH on heart function, but there are 
some limitations to our work. Our ability to investigate the growth and remodelling 
was limited by the fact that we didn’t account for the three-dimensional structural 
effects from the heart chambers, which alter the anatomy and behavior of the RV. To 
overcome this limitation, our framework could incorporate a 3D finite element model 
of the ventricles [21]. This extension of the lumped-parameter model would provide 
a more detailed characterization of ventricular mechanics as well as a more realistic 
representation of ventricular interaction via the septum. On the other hand, the limited 
literature sources with specific experimental data relevant to our study partly affect 
the accuracy of our results. There is significant variability in the biomechanical 
data between different cardiac studies, which makes it challenging to incorporate 
elements from multiple studies and unify them in one model [8, 11]. Conducting our 
own in vivo experiments could be useful to obtain more reliable parameters for the 
model. That is, experimental measurements would be a key next step to corroborate 
the results and develop an improved version of our model. 


7.5 Conclusion 


In this computational study, we have implemented a closed-loop lumped-parameter 
model of the heart and circulatory system for studying the cardiovascular adaptations 
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in the RV to PAH. We demonstrated that RVAD and atrial septostomy are effective 
treatments to improve left ventricular filling and cardiac output. We also extended 
our base model to account for direct ventricular interaction via the interventricular 
septum and the presence of the pericardium, showing that the pericardium has 
a larger effect relative to the ventricular interaction. Our present implementation 
lays the foundation for the addition of more components to the model that would 
augment its realism even further. We believe that highly complex computational 
models will continue to be critical to achieve our vision of dynamically simulating 
and characterizing ventricular mechanics in PAH. 
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